向量化循环答案

【问题标题】：Vectorizing a loop向量化循环
【发布时间】：2011-07-03 22:39:07
【问题描述】：

我正在创建一些人工数据。我需要创建家庭ID（H_ID）和个人ID（P_ID，在每个家庭）。

我找到了一种以矢量化方式创建 H_ID 的方法。

N <- 50

### Household ID
# loop-for
set.seed(20110224)
H_ID <- vector("integer", N)
H_ID[1] <- 1
for (i in 2:N) if (runif(1) < .5) H_ID[i] <- H_ID[i-1]+1 else H_ID[i] <- H_ID[i-1]
print(H_ID)

# vectorised form
set.seed(20110224)
r <- c(0, runif(N-1))
H_ID <- cumsum(r < .5)
print(H_ID)

但我不知道如何以矢量化方式创建 P_ID。

### Person ID
# loop-for
P_ID <- vector("integer", N)
P_ID[1] <- 1
for (i in 2:N) if (H_ID[i] > H_ID[i-1]) P_ID[i] <- 1 else P_ID[i] <- P_ID[i-1]+1
print(cbind(H_ID, P_ID))

# vectorised form
# ???

【问题讨论】：

标签： r loops for-loop

【解决方案1】：

另一个例子：

P_ID <- ave(rep(1, N), H_ID, FUN=cumsum)

几天前我发现了ave 函数（此处），并发现它在许多情况下都是非常有用且高效的快捷方式。

【讨论】：

【解决方案2】：

P_ID <- unname(unlist(tapply(H_ID, H_ID, function(x)c(1:length(x)))))

【讨论】：

【解决方案3】：

受Martin Morgan's solution 对一个密切相关的问题的启发，这是使用cummax 函数生成P_ID 的真正矢量化方法。一旦您注意到P_ID 与!(r < 0.5) 的cumsum 密切相关，就会很清楚：

set.seed(1)
N <- 10
r <- c(0, runif(N-1))
H_ID <- cumsum(r < .5)
r_ <- r >= .5 # flip the coins that generated H_ID.
z <- cumsum(r_)  # this is almost P_ID; just need to subtract the right amount...
# ... and the right amount to subtract is obtained via cummax
P_ID <- 1 + z - cummax( z * (!r_) )
> cbind(H_ID, P_ID)
      H_ID P_ID
 [1,]    1    1
 [2,]    1    2
 [3,]    2    1
 [4,]    3    1
 [5,]    3    2
 [6,]    3    3
 [7,]    3    4
 [8,]    4    1
 [9,]    5    1
[10,]    5    2

我没有做过详细的时序测试，但它可能很快就坏了，因为这些都是内部的矢量化函数

【讨论】：

我进行了计时测试 (N <- 2e6)。您的解决方案肯定是最快的。与lapply 解决方案相比，它快了大约 34 倍。谢谢！

【解决方案4】：

seq_along() 在这里是一个有用的工具。此示例将 H_ID 自身拆分为包含家庭的列表：

> head(split(H_ID, H_ID))
$`1`
[1] 1 1

$`2`
[1] 2

$`3`
[1] 3 3 3 3
....

然后，Q 的解决方案是将lapply() seq_along() 函数用于每个列表元素； seq_along() 创建一个向量 1:length(foo)。最后两个管家步骤，取消列出结果然后删除names：

> unname(unlist(lapply(split(H_ID, H_ID), seq_along)))
 [1] 1 2 1 1 2 3 4 1 1 2 3 1 1 1 1 1 2 3 4 5 1 2 3 4 1 1 2 1 2 1
[31] 1 2 1 2 3 4 1 2 1 2 1 2 1 1 2 1 2 1 2 3

【讨论】：

【解决方案5】：

这是一个相当紧凑且富有表现力的解决方案。就其中间值而言，与辛普森的有些相似：

cbind(H_ID,   unlist( sapply(table(H_ID), seq) ) )

其策略的核心是使用 table()-ed 值作为 seq() 的参数，默认情况下将采用单个数值并从 1 返回一个序列。

【讨论】：