ベクトル化されループ化されたバージョンが異なる回答を返す

Question

私は製品拡散シミュレーション研究を行っています。シミュレーションはノードのネットワークから始まり、小さな初期数のノードを製品にシードします。シーディング段階を超えた拡散は、製品を採用したノードの隣接ノードの数に依存する確率ルールによって管理されます。私はRでこのモデルの2つのバージョンを書いています - 1つはベクトル化され、もう1つはループ付きです。このアイデアはおそらくコードで表現されています。このコードを実行

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ループと「ベクター化」バージョンでは、2つの非常に異なる処理が行われています。 forループの1000回の反復のすべてで、活性化確率ベクトルを更新します（常に増加しているか、少なくとも減少していません）。ベクトル化されたバージョンでは、1000回の反復がすべて「同じ時間に」実行されるので、確率はすべて同じベクトルを使用して計算されます。

たとえば、forループの最初の反復で、最初のノードが更新される確率を計算します。そうであれば、それは第2ノードが更新される確率に影響を及ぼす可能性がある。ベクトル化されたバージョンでは、第2のノードが更新される確率は、第1のノードが更新されるかどうかに影響されない。

両方を同じにしたい場合は、新しいステータスを計算するときに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を取得します。しかし、同じことをするためには、実行する前に種を設定する必要があります。