代码之家 › 专栏 › 技术社区 › RMacey

R中的矢量化(子集)赋值

vectorization subset r

RMacey · 技术社区 · 8 年前

我对以下代码的结果感到惊讶。我希望(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将由sample函数生成(替换是索引可能重复的原因),w的长度将更长。这将是蒙特卡罗模拟多次运行的一部分,因此我希望它速度更快,从而避免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 比使用table函数快得多。

2 回复 | 直到 8 年前

Martin Morgan 8 年前

这会很快

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

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

以下是作为函数的解决方案

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

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

r2evans 8 年前

您看到的这种索引行为通常是需要的,特别是在类似“字典查找”的场景中,您希望查找一次,然后将其与其他场景分开保存。这是穷人的“加入”或“合并”操作:

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

不幸的是,您需要找到一种方法来迭代每个值,并以其他方式完成其工作。

人们可能会想试试 sapply 或者它的一个朋友,但一个计算的状态不执行:每次函数(的第二个参数 多愁善感的 )调用时,它不知道上次返回的内容。

所以你需要做一个类似的滚动动作。您可以使用 zoo::rollapply ,但另一种技术是“减少”,其中前一步的返回值是此迭代的输入。我们将初始条件设置为原始零向量 w ,并在每个 subw :

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

这实际上是在呼叫

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

出于教学目的,您可以将其汇总到:

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

(顺便说一句:引擎盖下, Reduce 实际上是使用 for 循环,但我更喜欢使用它,因为它清楚(至少对我来说)正在发生什么。加码高尔夫。)