Python实现新型冠状病毒传播模型及预测代码实例

1.传染及发病过程

一个健康人感染病毒后进入潜伏期（时间长度为Q天），潜伏期之后进入发病期（时间长度为D天），发病期之后该患者有三个可能去向，分别是自愈、接收隔离、死亡。

2.模型假设

潜伏期Q=7天，根据报道潜伏期为2~14天，取中间值；发病期D=10天，根据文献报告，WHO认定SARS发病期为10天，假设武汉肺炎与此相同；潜伏期的患者不具有将病毒传染给他人的能力；发病期的患者具有将病毒传染给他人的能力；患者在发病期之后不再具有将病毒传染他人的能力；假设处于发病期的患者平均每天密切接触1人，致使该人患病的概率为γ最初只有一个人类感染者；病情自然发展，没有外部干扰。

3.模型公式

：人类感染该病毒的天数

: 第N天感染该病毒并且处于发病期的患者数量

： Q天前新被感染患者，当日进入发病期的数量

： 当日发病期满，不再具有传染能力的患者数量

4.模型初始值

根据假设，最初只有一个人类感染者，所以：

5.实际疫情数据

人类感染病毒且发病的初始日期：根据财新网的报道，官方通报首例不明原因肺炎是在12月8日，考虑到确诊之前肯定已经尝试过各种治疗方案无效后认定为不明原因肺炎，所以有理由认为该名患者在12月8日已经处于发病期末端，根据假设发病期为10天，所以可以假设该名患者在11月29号发病，即N=1对应11月29日。

日新增发病数

近期疫情防控大事记：

从以上信息可以判断核酸检测试剂是在1月16日、17日大幅使用的，18日、19日确诊大量病例，因此1月20日之前的确诊病例数对模型参考意义不大。1月20日之后，可以认为新发病例，发病即检测。

截至1月20日24时，国家卫健委公告累计确诊病例291

截至1月21日24时，国家卫健委公告累计确诊病例440

截至1月22日24时，国家卫健委公告累计确诊病例571

截至1月23日24时，国家卫健委公告累计确诊病例830

截至1月24日24时，国家卫健委公告累计确诊病例1287

所以：

1月21日新增确诊病例：440-291=1491月22日新增确诊病例：571-440=1311月23日新增确诊病例：830-571=2591月24日新增确诊病例：1287-830=457

考虑到：

1月20日之后发病即检测确诊检测用时2日国外新增病例在个位数，且不能保证发病即检测等因素，暂不考虑新增病例应该递增，所以1月22日新增数据异常，舍去

使用上溯每日新增数据,同时考虑到确诊需要2天， 可以得到:

, ,

（ 对应11月29日 ）

6.拟合确定

根据近期每日新增数据、模型初始值及模型公式，用最小二乘拟合得到

7.预测

患者数量：根据上文确定的模型及参数，从11月29日（N=1）至1月27日（N=60）人群中累计处于发病期的人数如下图所示：

根据模型，近期人群中患者数量计算如下：

人群中感染了病毒并处于发病期的患者数量，注意Pn一般需要延期2日才能确诊

每日新增患者数量：

根据模型，近期人群中每日新增患者数量计算如下：

注意：图中是人群中新增发病患者数量，可与晚2日的政府发布新增数量进行对比。即22日新增患者数量可与24日政府发布的新增病例进行对比。截至目前模型计算22日新增为369人，政府公布的24日新增病例457人

根据以上模型预计未来几日的情况如下：

1月25日将新确诊432例，人群中发病患者为4068人；

1月26日将新确诊505例，人群中发病患者为4759人；

1月27日将新确诊590例，人群中发病患者为5568人；

1月28日将新确诊691例，人群中发病患者为6514人；

1月29日将新确诊809例，人群中发病患者为7621人。

由于1月20日之后采取的各种措施，将导致发病期D下降，感染概率 下降，1月29日之后日均增长势头会减弱。

8.模型代码

import numpy as np
import matplotlib.pyplot as plt
 
gamma = 0.55
Q = 7
D = 10
P = np.zeros(300, dtype=np.float)
Psum = np.zeros(300, dtype=np.float)
for i in range(Q):
  P[i] = 1
for j in range(300-Q):
  P[j+Q] = P[j+Q-1]+P[j]*gamma
  if j+Q-D-Q >= 0:
    P[j+Q] -= P[j+Q-D-Q]*gamma
  if j+Q == D:
    P[j+Q] -= 1
plt.xlabel("N")
plt.ylabel("PN")
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.plot(range(1, 61), P[0:60])
plt.grid()
plt.show()