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

所属分类: 脚本专栏 / python 阅读数: 2012
收藏 0 赞 0 分享

1.传染及发病过程

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

2.模型假设

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

3.模型公式

P_{N} = P_{N-1} + P_{N-Q} \cdot \gamma + P_{N-Q-D} \cdot \gamma

N:人类感染该病毒的天数

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

P_{N-Q} \cdot \gamma : Q天前新被感染患者,当日进入发病期的数量

P_{N-Q-D} \cdot \gamma: 当日发病期满,不再具有传染能力的患者数量

4.模型初始值

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

P_{1} = P_{2} = P_{3} ... = P_{7} = 1

5.实际疫情数据

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

日新增发病数 \Delta P_{n}

近期疫情防控大事记:

从以上信息可以判断核酸检测试剂是在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天, 可以得到:

\Delta P_{52} = 149, \Delta P_{54} = 259, \Delta P_{55} = 457

P_{1} 对应11月29日 )

6.拟合确定 \gamma

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

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下降,感染概率 \gamma 下降,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()

以上就是脚本之家小编给大家整理的全部内容,如果大家有任何补充可以联系小编。

更多精彩内容其他人还在看

pandas的qcut()方法详解

这篇文章主要介绍了pandas的qcut()方法详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

从列表或字典创建Pandas的DataFrame对象的方法

这篇文章主要介绍了从列表或字典创建Pandas的DataFrame对象的方法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

pandas.DataFrame的pivot()和unstack()实现行转列

这篇文章主要介绍了pandas.DataFrame的pivot()和unstack()实现行转列,小编觉得挺不错的,现在分享给大家,也给大家做个参考。一起跟随小编过来看看吧
收藏 0 赞 0 分享

pandas中的series数据类型详解

这篇文章主要介绍了pandas中的series数据类型详解,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

pandas 时间格式转换的实现

这篇文章主要介绍了pandas 时间格式转换的实现,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

python中时间、日期、时间戳的转换的实现方法

这篇文章主要介绍了python中时间、日期、时间戳的转换的实现方法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

pandas进行时间数据的转换和计算时间差并提取年月日

这篇文章主要介绍了pandas进行时间数据的转换和计算时间差并提取年月日,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

详解将Pandas中的DataFrame类型转换成Numpy中array类型的三种方法

这篇文章主要介绍了详解将Pandas中的DataFrame类型转换成Numpy中array类型的三种方法,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧
收藏 0 赞 0 分享

python和c语言的主要区别总结

在本篇文章里小编给各位整理了关于python和c语言的主要区别的相关知识帖内容,有需要的朋友们学习阅读下。
收藏 0 赞 0 分享

选择Python写网络爬虫的优势和理由

在本篇文章里小编给各位整理了一篇关于选择Python写网络爬虫的优势和理由以及相关代码实例,有兴趣的朋友们阅读下吧。
收藏 0 赞 0 分享
查看更多