这个模拟很容易;不过我想再次证实两点关键的模拟规则:
(1)每次死亡的物种是随机选择的,与它们当前的物种数目无关;
(2)每次替代的物种是基于当前的物种数目的(以该数目为概率选取物种);
如果这两点没有问题的话,那么我估计5行左右的代码就可以完成模拟了,这里是一个图示:
ecol.simu = function(nr = 10, nc = 10, col.sp = c(1, <br />
2), pch.sp = c(1, 2), col.die = 1, pch.die = 4, cex = 3, <br />
nmax = 50, interval = 1) {<br />
x = rep(1:nc, nr)<br />
y = rep(1:nr, each = nc)<br />
par(ann = FALSE)<br />
p = sample(rep(1:2, nr * nc/2), nr * nc)<br />
for (i in 1:nmax) {<br />
plot(1:nc, 1:nr, type = "n", xlim = c(0.5, nc + 0.5), <br />
ylim = c(0.5, nr + 0.5))<br />
abline(h = 1:nr, v = 1:nc, col = "lightgray", lty = 3)<br />
points(x, y, col = col.sp[p], pch = pch.sp[p], cex = cex)<br />
Sys.sleep(interval)<br />
idx = sample(nr * nc, 1)<br />
points(x[idx], y[idx], pch = pch.die, col = col.die, <br />
cex = cex, lwd = 3)<br />
tbl = as.vector(table(p))<br />
tbl = tbl + sign(p[idx] - 1.5) * c(1, -1)<br />
p[idx] = sample(1:2, 1, prob = tbl)<br />
Sys.sleep(interval)<br />
}<br />
p<br />
}<br />
par(mar = c(3, 3, 1, 1))<br />
# 正常的模拟<br />
ecol.simu()<br />
# 让死亡来得更猛烈些吧!!<br />
ecol.simu(col.sp = c(8, 2), pch.sp = c(20, 17), nmax = 1000, <br />
interval = 0.05)