矢量化和循环版本返回不同的答案

Question

我正在运行产品扩散模拟研究。仿真从一个节点网络开始，并将一个初始节点数量与一个产品进行种子初始化。播种阶段之后的扩散受概率规则控制，该规则取决于采用产品的节点的邻居数量。我在R中编写了这个模型的两个版本 - 一个是矢量化的，另一个是循环的。这个想法可能更好地用代码表达。

library(igraph) 

set.seed(20130810) 

g <- sample_smallworld(dim = 1, size = 1000, nei = 12, p = 0.6) 

n.nodes <- length(V(g)) 
nbr.influence <- rnorm(n = n.nodes, mean = 0.18, sd = 0.01) 

# Diffusion simulation with loops 

nodes.status <- rep.int(0, n.nodes) 
seed <- sample(V(g), size = as.integer(0.005*n.nodes)) 
nodes.status[seed] <- 1 

cat("Number of nodes seeded (loop version): ", sum(nodes.status), "\n") 

for (node in V(g)) { 
    if (nodes.status[node] != 1) { 
    n.active.nbrs = 0 

    for (nbr in neighbors(g, node)) { 
     if (nodes.status[nbr] == 1) n.active.nbrs <- n.active.nbrs + 1 
    } 

    prob.change <- 1 - (1 - nbr.influence[node])^n.active.nbrs 

    if (runif(n = 1) < prob.change) nodes.status[node] = 1 
    } 
} 

cat("Number of nodes engaged after one iteration (loop version): ", 
sum(nodes.status), "\n") 

# Vectorized diffusion simulation 

A <- get.adjacency(g) 

nodes.status <- rep.int(0, n.nodes) 
seed <- sample(V(g), size = as.integer(0.005*n.nodes)) 
nodes.status[seed] <- 1 
cat("Number of nodes seeded (vectorized version): ", sum(nodes.status), "\n") 

# use the adjacency matrix to count number of active neighbors for each node 
n.active.neighbours <- as.vector(A %*% nodes.status) 

# build the activation probability vector 
prob.change <- 1 - (1 - nbr.influence)^n.active.neighbours 

# see which of the nodes are ready to activate 
vuln.nodes <- runif(n = n.nodes) < prob.change 

# activate those nodes which are ready 
nodes.status[vuln.nodes > nodes.status] <- 1 

cat("Number of nodes engaged after one iteration (vectorized version): ", 
sum(nodes.status), "\n") 

运行该代码给出以下输出

Number of nodes seeded (loop version): 5 
Number of nodes engaged after one iteration (loop version): 380 
Number of nodes seeded (vectorized version): 5 
Number of nodes engaged after one iteration (vectorized version): 32 

两个版本的逻辑是相同的（即，扩散是由相同的概率规则），但最终的答案是广泛不同。这段代码中的错误在哪里？

Answer 1

您的for循环和“矢量化”版本正在做两件完全不同的事情。在for循环的所有1000次迭代中，更新（始终增加或至少不减少）激活概率向量。在向量化版本中，所有1000次迭代都是“同时”完成的，所以概率都是使用相同的向量计算的。

例如，在for循环的第一次迭代中，计算第一个节点更新的概率。如果是这样，那么它会影响第二个节点更新的可能性。在矢量化版本中，不管第一个节点是否更新，第二个节点更新的概率都不受影响。

如果要使两者相同，则必须在计算新状态时使循环保持所有其他节点的原始状态。这里是一个例子：

cat("Number of nodes seeded (loop version): ", sum(nodes.status), "\n") 

new.nodes.status <- nodes.status # copy vector to preserve original state. 
for (node in V(g)) { 
    if (nodes.status[node] != 1) { 
    n.active.nbrs = 0 
    for (nbr in neighbors(g, node)) { 
     if (nodes.status[nbr] == 1) n.active.nbrs <- n.active.nbrs + 1 
    } 
    prob.change <- 1 - (1 - nbr.influence[node])^n.active.nbrs 
    if (runif(n = 1) < prob.change) new.nodes.status[node] = 1 # Only update new. 
    } 
} 

cat("Number of nodes engaged after one iteration (loop version): ", 
    sum(new.nodes.status), "\n") 

哪个让我32。但是你应该在每次运行之前设置你的种子以获得完全相同的结果。