Python基于SIR模型模拟预测新冠疫情数据

Theo

发布时间：2026-10-13 | 4人浏览过

原创

本文使用Python构建SIR传染病模型，通过手工提取数据并采用穷举法估算感染系数β和初始易感人数S0，结合已知恢复率γ进行疫情走势模拟。结果显示最大感染人数约2.3万人，峰值出现在2月12日左右。文章强调此为个人技术探索，非专业流行病学预测。

大家好。

SIR模型原理与假设

背景无需赘述，原定初四复工，现延期至2月10日。这是我职业生涯中最长的假期，无奈无法外出，便利用Python对新冠病毒肺炎数据进行了模拟预测。需声明，本文纯属个人兴趣研究，不代表客观事实。

采用SIR模型，其中S代表易感人群，I代表感染者，R代表恢复者。感染者作为传染源，以一定概率将病毒传给易感者，同时自身也有概率被治愈（获得免疫）或死亡。易感者一旦感染即转化为新的传染源。

模型假设如下：

第一，忽略人口出生、死亡及流动，假设总人口数量恒定。

第二，假设患者与易感者接触即具备传染性。设t时刻单位时间内，单个患者传染的易感者数量与该环境内易感者总数s(t)成正比，比例系数为β。因此，t时刻单位时间内所有患者传染的总人数为βs(t)i(t)。

第三，设t时刻单位时间内从感染者中移出（治愈或死亡）的人数与患者数量成正比，比例系数为γ，移出量为γi(t)。模型公式为：

其中，β为感染系数，表示易感者与传染源接触后被感染的概率；γ为隔离（恢复）系数，其倒数1/γ代表平均感染周期。S(0)为初始易感人数，I(0)为初始感染人数。

参考[1]中的代码，该模型下的感染人数变化如下：

数据获取与参数估算

当前关键在于利用现有数据确定新冠肺炎的β值和γ值。首先获取数据，从[3]中提取，因数据量较小，采用手工方式处理并保存为csv文件。

随后对数据进行可视化绘图。

另一个关键指标是基本再生数R0=β/γ，当R0大于1时，人群中大部分可能被感染[4]。世卫组织1月23日估计R0在1.4至2.5之间[5]，基于前425例数据的最新估计值为2.2[6]。

文献[7]按一般病毒性肺炎恢复期25天计算，得出γ值为0.04。

Python实现与预测结果

关于β值和初始易感人群，文献[7]采用先估计区间再用最小二乘法寻找最佳参数的方法，得出β≈3.57*10^-5，S[0]范围在5000至30000人之间。[7]提供了Matlab代码，我将其改写为Python。由于对最小二乘法实现不熟悉，经过多次尝试，最终决定采用最直接的穷举法。即通过两层嵌套循环遍历范围内所有β值和S0值，计算每次尝试结果与实际数据差值的平方和，选取平方和最小的一组参数进行预测。代码如下：

设定γ值为0.04，对应一般病程25天。

利用穷举法估算β值和初始易感人数。


gamma = 0.04
S0 = [i for i in range(20000, 40000, 1000)]
beta = [f for f in np.arange(1e-7, 1e-4, 1e-7)]
# 定义偏差函数
def error(res):
 err = (data["感染者"] - res)**2
 errsum = sum(err)
 return errsum

# 穷举法，找出与实际数据差的平方和最小的S0和beta值
minSum = 1e10
minS0 = 0.0
minBeta = 0.0
bestRes = None

for S in S0:
 for b in beta:
  # 模型的差分方程
  def diff_eqs_2(INP, t):
   Y = np.zeros((3))
   V = INP
   Y[0] = -b * V[0] * V[1]
   Y[1] = b * V[0] * V[1] - gamma * V[1]
   Y[2] = gamma * V[1]
   return Y

  # 数值解模型方程
  INPUT = [S, I0, 0.0]
  RES = spi.odeint(diff_eqs_2, INPUT, t_range)
  errsum = error(RES[:21, 1])
  if errsum < minSum:
   minSum = errsum
   minS0 = S
   minBeta = b
   bestRes = RES
   print("S0=%d beta=%f minErr=%f" % (S, b, errsum))
print("S0 = %d β = %f" % (minS0, minBeta))

计算结果为 S0 = 39000, β = 8e-6。

上述程序运行耗时较长，仅在探索阶段执行，确定最优参数后将其注释，改用最优参数进行正式预测。