R语言：向量计算中的随机赋值陷阱

在 R 语言的数据模拟中，向量计算往往因其高效性而被广泛采用。当我们使用向量计算来模拟个体状态的演变时，如果不加注意，可能会无意间引入模式化的行为，影响模拟的随机性和异质性。

向量计算的吸引力与隐患

在 R 语言中，向量化操作比 for 循环更快，这是因为 R 的底层计算基于 C 语言优化，使得整个向量可以一次性传递给计算函数，而无需逐个元素进行解释和处理。

在根据规则来生成模拟数据时，通常优先考虑向量计算（ifelse），而非循环遍历。然而，在涉及随机赋值的情境下，向量计算可能会意外地引入“锁步”的陷阱。

同步演变下的异常模式

考虑这样一个模拟场景，每个个体的状态数值可以在一定范围内波动，遵循以下规则：

若状态为正数，则逐步减少；若状态为负数，则逐步增加。
若状态为 0，则该个体会被赋予一个新的随机值。

ifelse(state == 0, sample(new_values, 1), state ± 1)

这个过程应该保持较高的个体异质性。当我们使用向量计算直接实现时，个体状态是同步变化。

set.seed(1)

ca <- data.frame('step0' = round(runif(50, min = -100, max = +100)))

for (t in 1:100) {
  ca_rule <- ifelse(ca[, t] == 0, sample(-100:100, 1),
                    ifelse(ca[, t] > 0, ca[, t] - 1, ca[, t] + 1))
  ca[,ncol(ca) + 1] <- ca_rule
  colnames(ca)[ncol(ca)] <- paste0("step", t)
}

rownames(ca) <- paste0("Student ", 1:50)

表面上看，这段代码完美执行了我们想要的规则。然而，这里的问题在于 sample(-100:100, 1) 仅生成一个随机数，然后将这个单一数值赋给所有状态为 0 的个体。

如果某个时间步内有多个个体同时满足 state == 0 的条件，它们将在该时间步获得完全相同的值，并在接下来的演化中保持相同的变化趋势。

pheatmap(ca, 
         cluster_rows = TRUE, 
         cluster_cols = FALSE,
         show_rownames = TRUE,
         show_colnames = FALSE,
         treeheight_row = 0)

这种模式的典型特征是：

如上图，某些个体状态的演化轨迹变得完全相同，形成“锁步”。
预期的个体异质性被破坏，可能导致整个模拟结果失真。
若同步个体的比例较大，模拟结果可能与预期目标产生较大偏差。

这种现象削弱了模拟的随机性，影响个体层面的异质性，甚至改变原本预期的模拟行为。

如何避免？

要保持个体间的独立性，我们需要确保每个个体的状态更新过程真正独立，特别是在随机赋值时。最直接的方法是为每个需要重新赋值的个体单独生成随机数，而不是使用一个统一的随机值。

一种方法是先找到所有满足 state == 0 的个体索引，然后逐个赋予不同的随机值。这一方式可以避免同步陷阱，同时仍然保持向量计算的高效性。

vectorized_update <- function() {
  ca <- data.frame('step0' = round(runif(50, min = -100, max = +100)))

  for (t in 1:1000) {
    zero_indices <- which(ca[, t] == 0)  
    new_random_values <- sample(-100:100, length(zero_indices), replace = TRUE)  

    ca_rule <- ifelse(ca[, t] == 0, NA,  
                      ifelse(ca[, t] > 0, ca[, t] - 1, ca[, t] + 1))
    
    if (length(zero_indices) > 0) {
      ca_rule[zero_indices] <- new_random_values
    }
    
    ca[, ncol(ca) + 1] <- ca_rule
    colnames(ca)[ncol(ca)] <- paste0("step", t)
  }
}

与之相比，传统的循环遍历方式可以逐个为每个个体生成随机数，从而保证了个体更新的独立性，避免了“锁步”模式的出现，但这种方法往往牺牲了运行效率。

loop_update <- function() {
  ca <- data.frame('step0' = round(runif(50, min = -100, max = +100)))

  for (t in 1:1000) {
    ca_rule <- numeric(50)
    for (i in 1:50) {
      if (ca[i, t] == 0) {
        ca_rule[i] <- sample(-100:100, 1)
      } else if (ca[i, t] > 0) {
        ca_rule[i] <- ca[i, t] - 1
      } else {
        ca_rule[i] <- ca[i, t] + 1
      }
    }
    ca[, ncol(ca) + 1] <- ca_rule
    colnames(ca)[ncol(ca)] <- paste0("step", t)
  }
}

ifelse() 向量化操作的底层用 C 语言实现，比 for 逐个遍历快很多。
同时，因为 sample() 在 vectorized_update() 中一次性生成多个随机数，而 loop_update() 逐个生成，调用次数远多于向量化方案。
此外，which(ca[, t] == 0) 一次性找到所有 state == 0，也比循环判断更高效。

预期向量化方案将比循环遍历快得多，尤其当 t 变大时，累积性能差距更明显。

replicate(5, system.time(vectorized_update())["elapsed"])
0.340000000000146
0.319999999999709
0.340000000000146
0.270000000000437
0.25

replicate(5, system.time(loop_update())["elapsed"])
2.5
2.42000000000189
2.40999999999985
2.90000000000146
2.88999999999942

通过 system.time() 测量后，我们可以发现向量化方法在大规模数据下通常能显著缩短运行时间，而循环遍历的开销较大。

结语

尽管向量化计算在性能上具有无可争议的优势，尤其在大规模数据或长时间步模拟中能够显著缩短运行时间，但在处理随机化操作时必须格外谨慎。

对比之下，采用逐个遍历的方式（即使用 for 循环逐个处理每个元素）虽然在性能上有所牺牲，却能确保每个元素获得独立的随机数，从而保留向量中各个元素的独立性，进而提升模拟结果的科学性和可靠性。

#R #Data Visualization #Complexity #Agent-Based Model ...all tags