netlogo - SIR 模型错误 - 找不到错误,需要帮助找出潜在的偏差源?
问题描述
这个问题将是一个有趣的问题。我试图在一个自由移动的代理系统中复制一篇关于疾病传播的论文的结果(听起来像是 NetLogo 的完美工作)。根据论文中给出的细节,我很容易在 NetLogo 中编写了一个简单的 SIR 模型,确保我的模型参数与列出的参数相匹配,然后让模拟运行。一切都运行得很完美,直到我检查了实验结果与预测值的匹配情况(根据论文的结果)。他们离开了,而且差距相当大。认为代码中的某个地方有错误,我检查了所有内容,但什么也没找到. 然后我确保事件的顺序是正确的(因为移动、感染和恢复的顺序很重要),这些也与论文相符。我仔细考虑了这个问题很长一段时间,直到最后我打开了 R,在 RStudio 中编写了完全相同的程序,然后让它运行,结果发现结果完全符合预测!R 代码和我期望NetLogo 代码做的事情一样,所以我认为 NetLogo 的幕后正在发生一些事情,或者我在某个地方有误解,这是偏差的根源......请注意,因为论文中的结果是平均场近似值,您必须运行该程序几次才能使其接近理论结果。
我不确定我哪里出错了,因为我的 R 代码确认预测值是正确的,所以我得出结论,我的 NetLogo 代码中的某个地方是不正确的。我对 NetLogo 不太熟悉,如果有人能帮助我找到以下代码中可能发生偏差的位置,我将不胜感激。实验平均值往往低于预测值,这表明感染发生的速度比预期的要快,但在我看到的所有变化中,没有一个能解决这个问题(例如,每只感染龟一次不会发生一次感染) . 任何建议/帮助将不胜感激。
下面介绍了我的代码的精简版本。这应该在带有标准设置/开始按钮的常规界面中运行。结果存储在可以绘制的列表中,任何好奇的人都可以通过 Plot 对象在模拟进行时看到偏差。先感谢您。
;; Simple SIR model
globals [
;; variables for storing predictions
predS
predE
predI
predR
oldPredS
oldPredE
oldPredI
oldPredR
;; list to store experimental values
Slist
;; list to store predicted values
predSList
;; model variables
length-of-patch ;; length of habitat (a square of area length-of-patch^2)
infection-radius ;; the distance from an infectious individual a susceptible agent has to be within
;; in order to risk getting infected
total-pop ;; total population in the model
force-of-infection ;; probability of infection if within infection-radius distance
I0 ;; initial infected
recovery-rate ;; probability of recovery
]
turtles-own [
infected-status ;; 0 susceptible, 1 infected, 2 recovered
]
to setup
ca ;; clear
;; define the variables
set length-of-patch 31.62278 ;; the square root of 1000 (so the density is 1)
set infection-radius 1
set total-pop 1000
set force-of-infection 0.1
set I0 10
set recovery-rate 0.05
;; setup simulation
setup-patches
setup-agents
reset-ticks
;; initialize lists as empty
set Slist []
set predSList []
end
to go
;; update experimental values (density of susceptible individuals)
set Slist lput ((count turtles with [infected-status = 0]) / (length-of-patch ^ 2)) Slist
if (ticks = 0) ;; if ticks == 0, make sure initial value is the same as experimental
[
;; update predicted values with densities of agents
set predS ((count turtles with [infected-status = 0]) / (length-of-patch ^ 2))
set predI ((count turtles with [infected-status = 1]) / (length-of-patch ^ 2))
set predR 0
;; placeholder variables for iterative process
set oldPredS predS
set oldPredI predI
set oldPredR predR
;; store predicted S population in corresponding list
set predSList lput (predS) predSList
]
if (ticks > 0) ;; if ticks > 0, then update predicted values according to paper results
[
;; update predicted values
set predI (oldPredI + oldPredS * (1 - (1 - force-of-infection * oldPredI) ^ (pi * (infection-radius ^ 2))) - recovery-rate * oldPredI)
set predR (oldPredR + recovery-rate * oldPredI)
set predS ((total-pop / (length-of-patch ^ 2)) - predI - predR)
;; placeholder variables
set oldPredS predS
set oldPredI predI
set oldPredR predR
;; store values in corresponding list
set predSList lput (oldPredS) predSList
]
;; perform movement, infection, and recovery, in that order
move-agents
infect-agents
recover-agents
if (count turtles with [infected-status = 1] = 0) [
;; if no one else is infected, stop
stop
]
tick
end
to setup-patches
;; resize the world to make it fit comfortably in the interface
resize-world 0 length-of-patch 0 length-of-patch
set-patch-size 400 / (length-of-patch)
end
to setup-agents
;; create susceptible agents
crt (total-pop - I0) [
set infected-status 0
setxy random-pxcor random-pycor
set color 55 ;; green
set size 2
]
;; create I0 infected agents
crt I0 [
set infected-status 1
setxy random-pxcor random-pycor
set color 15 ;; red
set size 2
]
end
to move-agents ;; move all the agents
ask turtles [
setxy random-pxcor random-pycor
]
end
to infect-agents
;; iterate over infected turtles
ask turtles with [infected-status = 1] [
;; check neighborhood around infected turtle for susceptible turtles...
let numNeighbors count (turtles with [infected-status = 0] in-radius infection-radius)
if (numNeighbors > 0) [ ;; there are susceptibles around, so we perform infection
ask (turtles with [infected-status = 0] in-radius infection-radius) [
let %draw (random-float 1)
if (%draw <= force-of-infection) [ ;; probability of infection
;; infect one of the neighbors
set infected-status 1
set color 15 ;; red
]
]
] ;; end of if numneighbors > 0
]
end
to recover-agents
ask turtles with [infected-status = 1] [
let %draw (random-float 1)
if (%draw <= recovery-rate) [ ;; an agent recovered
set infected-status 2
set color 105
]
]
end
解决方案
我可以看到的一个问题是你有:setxy random-pxcor random-pycor
但你想要:setxy random-xcor random-ycor
基本上,您将所有海龟放在补丁的中心,因此它们彼此重叠,而不是在空间中随机分布它们。这种定位改变了海龟之间可能距离的分布。
我还将海龟的数量更改为1024 1089,将大小更改为 sqrt 1024(而不是 1000),以使密度正确匹配。
两者都减少了不匹配,但由于我没有进行大量运行,因此尚不清楚它们是否解决了问题。
更新
甚至需要更多的维度匹配。更改代码以便有 1089 个代理,将长度设置为 33 以进行预测计算,并将世界大小调整为最大 32 似乎使曲线更接近。这认识到补丁坐标 0 到 32 实际上描述了长度为 33 的大小,因为 NetLogo 坐标将从 -0.5 开始并运行到 @Jasper 提到的 32.5
推荐阅读
- reactjs - 类型 'MainMenu[]' 不可分配给类型 'never[]'
- shell - 仅删除第一行和最后一行的尾随逗号
- android - 我有一个相机抖动项目,我在制作自定义预览相机时遇到了麻烦,这里有没有人可以帮助我应该使用什么包?
- caching - 替换时可能会丢失 Apollo Client Cache 数据
- javascript - 通过函数数组运行一个值
- java - 用 $1 替换正则表达式中的第一个 - JAVA
- azure - azure vm 的诊断设置
- java - Java - HttpClient 不释放线程
- python - python discord bot on_disconnect 消息
- javascript - 无法加载资源:net::ERR_FILE_NOT_FOUND,未捕获(承诺中)类型错误:无法获取