SIR模型错误-找不到错误，需要帮助来找到潜在的偏差源吗？

问题描述

这个问题将是一个有趣的问题。我试图复制paper的结果，该结果涉及在自由移动的代理系统中的疾病传播（听起来像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个代理，将pred计算的长度设置为33，并且将世界的大小调整为32（最大值）似乎会使曲线更接近。这可以识别出补丁坐标0到32实际上描述了一个长度为33的大小，因为NetLogo坐标将以-0.5开头并运行到32.5，如@Jasper所述