首页 > 解决方案 > 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

标签: netlogo

解决方案


我可以看到的一个问题是你有: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


推荐阅读