R矩阵值回收？答案

【问题标题】：R matrix values recycling?R矩阵值回收？
【发布时间】：2021-03-19 11:55:54
【问题描述】：

我正在关注这个关于 R https://www.youtube.com/watch?v=uAeSykgXnhg 中基于代理的建模的 YouTube 教程。

我使用自己的变量名复制了代码（这有助于我更好地理解导师的代码）。目的是跟踪人们在接触他人时如何感染covid-19。并非每次接触都会导致感染。在连续的模型运行中，感染人数=人口规模，未感染人数应为0。这是我复制的代码：

# define first agent
agents <- data.frame(agent_no = 1,
                     state = "e",
                     mixing = runif(1,0,1))

# specify agent population
pop_size <- 100

# fill agent data
for(i in 2:pop_size){
  agent <- data.frame(agent_no = i,
                      state = "s",
                      mixing = runif(1,0,1))
  agents <- rbind(agents, agent)
}

# specify number of model runs
n_times <- 10

# initialise output matrix 
out <- matrix(0, ncol = 2, nrow = n_times)

# run simple agent-based model
for(k in 1:n_times){
  for(i in 1:pop_size){
    # likelihood to meet others
    likelihood <- agents$mixing[i]
    # how many agents will they meet (integer). Add 1 to make sure everybody meets somebody
    connect_with <- round(likelihood * 3, 0) + 1 
    # which agents will they probably meet (list of agents)
    which_others <- sample(1:pop_size, 
                           connect_with, 
                           replace = T, 
                           prob = agents$mixing)
    for(j in 1:length(which_others)){
      contacts <- agents[which_others[j],]
      # if exposed, change state
      if(contacts$state == "e"){
        urand <- runif(1,0,1)
        # control probability of state change
        if(urand < 0.5){
          agents$state[i] <- "e"
        }
      }
    }
  }
  out[k,] <- table(agents$state)
}

查看输出时，一旦每个人都被感染（第一列），未感染人数（第二列）应该是 0，但我得到 100，我怀疑这是由于回收。

     [,1] [,2]
 [1,]   12   88
 [2,]   33   67
 [3,]   69   31
 [4,]   86   14
 [5,]   92    8
 [6,]   95    5
 [7,]   97    3
 [8,]   98    2
 [9,]   99    1
[10,]  100  100

我运行了一些诊断程序来查看发生了什么：

table(agents$state)
      e 
    100 

agents[agents$state == "s",]

    [1] agent_no state    mixing  
    <0 rows> (or 0-length row.names)

我认为 0-length row.names 是我的问题。结果应该是这样的：

     [,1] [,2]
 [1,]   12   88
 [2,]   33   67
 [3,]   69   31
 [4,]   86   14
 [5,]   92    8
 [6,]   95    5
 [7,]   97    3
 [8,]   98    2
 [9,]   99    1
[10,]  100    0

有人可以解释我做错了什么吗？非常感谢。

【问题讨论】：

为了使这个可重现，顶部需要一个 set.seed 语句。当我尝试使用 set.seed(123) 时，out[10, 2] 的值为 1。

标签： r matrix agent

【解决方案1】：

我将n_times 增加到 10000 并且找不到回收的证据。虽然这并不意味着它没有发生，但这意味着不幸的是，如果没有明确的设置，我们将无法重现该问题。所以我在这里的建议是未经证实的。

选项 1

鉴于您发现了一个以所有agents$state == "e" 结尾的场景，那么我将建议一个技巧，该技巧将始终找到至少一个"s"（实际上是您所知道的每个值之一）：

  out[k,] <- table(c("e", "s", agents$state)) - 1

我假设唯一可能的值是"e" 和"s"；如果还有其他的，这种技术完全依赖于我们确保每个可能的值至少被看到一次的前提，然后递减所有值。由于我们为每个可能的值“添加一个观察值”，因此从表中减去一个是安全的。有了这个技巧，你的支票应该是

table(agents$state)
#       e 
#     100 
table(c("e", "s", agents$state))
#       e       s 
#     101       1
table(c("e", "s", agents$state)) - 1
#       e       s 
#     100       0

因此回收不应该是一个因素。

选项 2

另一种更健壮的技术（即不需要包含所有可能的值）是强制长度，假设我们确定它应该是什么（我认为我们在这里做）：

z <- table(agents$state)
z
#   s 
# 100 
length(z) <- 2
z
#   s     
# 100  NA

由于您“知道”长度应始终为 2，因此您可以在其中硬编码 2。

选项 3

这个方法更加健壮一点，你不需要知道绝对长度，它们都会被扩展到最长返回的长度。

首先，可重现的样本数据：

set.seed(2021)
agents <- data.frame(agent_no = 1,
                     state = "e",
                     mixing = runif(1,0,1))
# specify agent population
pop_size <- 100
# fill agent data
for(i in 2:pop_size){
  agent <- data.frame(agent_no = i,
                      state = "s",
                      mixing = runif(1,0,1))
  agents <- rbind(agents, agent)
}
head(agents)
#   agent_no state    mixing
# 1        1     e 0.4512674
# 2        2     s 0.7837798
# 3        3     s 0.7096822
# 4        4     s 0.3817443
# 5        5     s 0.6363238
# 6        6     s 0.7013460

替换您的 for 循环：

for (k in 1:n_times) {
}

与

out <- lapply(seq_len(n_times), function(k) {
  for(i in 1:pop_size){
    # likelihood to meet others
    likelihood <- agents$mixing[i]
    # how many agents will they meet (integer). Add 1 to make sure everybody meets somebody
    connect_with <- round(likelihood * 3, 0) + 1 
    # which agents will they probably meet (list of agents)
    which_others <- sample(1:pop_size, 
                           connect_with, 
                           replace = T, 
                           prob = agents$mixing)
    for(j in 1:length(which_others)){
      contacts <- agents[which_others[j],]
      # if exposed, change state
      if(contacts$state == "e"){
        urand <- runif(1,0,1)
        # control probability of state change
        if(urand < 0.5){
          agents$state[i] <- "e"
        }
      }
    }
  }
  table(agents$state)
})

此时，您有一个列表，可能是长度为 2 的向量：

out[1:3]
# [[1]]
#  e  s 
#  1 99 
# [[2]]
#  e  s 
#  2 98 
# [[3]]
#  e  s 
#  3 97

请注意，我们可以确定所有这些的长度

lengths(out)
#  [1] 2 2 2 2 2 2 2 2 2 2

类似于我们强制向量长度的选项 2，我们可以在这里做同样的事情：

maxlen <- max(lengths(out))
out <- lapply(out, `length<-`, maxlen)
## or more verbosely
out <- lapply(out, function(vec) { length(vec) <- maxlen; vec; })

您可以用table(lengths(out)) 确认它们的长度相同，应该是2 by n_times of 10。

从这里，我们可以将所有这些向量组合成一个矩阵

out <- do.call(rbind, out)
out
#        e  s
#  [1,]  1 99
#  [2,]  2 98
#  [3,]  3 97
#  [4,]  2 98
#  [5,]  1 99
#  [6,] 20 80
#  [7,] 12 88
#  [8,]  1 99
#  [9,]  2 98
# [10,]  1 99

【讨论】：