python模拟预测一下新型冠状病毒肺炎的数据
作者:自由民 发布时间:2023-12-17 05:09:15
大家还好吗?
背景就不用多说了吧?本来我是初四上班的,现在延长到2月10日了。这是我工作以来时间最长的一个假期了。可惜哪也去不了。待在家里,没啥事,就用python模拟预测一下新冠病毒肺炎的数据吧。要声明的是本文纯属个人自娱自乐,不代表真实情况。
采用SIR模型,S代表易感者,I表示感染者,R表示恢复者。染病人群为传染源,通过一定几率把传染病传给易感人群,ta自己也有一定的几率被治愈并免疫,或死亡。易感人群一旦感染即成为新的传染源。
模型假设:
①不考虑人口出生、死亡、流动等情况,即人口数量保持常数。
②一个病人一旦与易感者接触就必然具有一定的传染力。假设 t 时刻单位时间内,一个病人能传染的易感者数目与此环境内易感者总数s(t)成正比,比例系数为β,从而在t时刻单位时间内被所有病人传染的人数为βs(t)i(t)。
③ t 时刻,单位时间内从染病者中移出的人数与病人数量成正比,比例系数为γ,单位时间内移出者的数量为γi(t)。
模型为
其中,β为感染系数,代表易感人群与传染源接触被感染的概率。γ为隔离(恢复)系数,我们对其倒数1/γ更感兴趣,代表了平均感染时间(average infectious period)。S(0)为初始易感人数,I(0)为初始感染人数。
按照[1]里面的代码模型的感染人数是这样的
现在的问题就是利用现有的数据找到新冠肺炎的β值,γ值等数据了。先把数据拔下来吧。从[3]上扒数据,由于数据不多,就手工完成吧。保存到csv文件里。
然后把数据作图
还有一个指标是再生数R0=β/γ,大于1时人群中大部分才被感染[4]。世卫组织1月23日的估计是R0在1.4到2.5之间[5],最新的根据前425例发病数据的估计值为2.2[6]。
文章[7]中的按一般病毒性肺炎恢复期25天计算得到的γ值为0.04。
关于β值和初始易感人群,[7]的作者采用的方法是先估计一个区间,然后用最小二乘法找到最佳参数,β≈3.57*10^-5。S[0]的范围为5000-30000人。[7]文章里有matlab代码,我用python改写一下,由于对最小二乘法法的实现比较陌生,尝试了半天,最后我决定用最笨的办法——穷举法。就是用两个嵌套循环将范围内所有β值和S0值都试一遍,计算每次尝试结果与实际数据之间差值的平方和,平方和最小的一组β值和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
上述程序耗时较长,只在探索时执行,完了就注释掉,用最优参数进行预测。
预测最大感染人数:23769 时间是在1月10日的33天后,也就是2月12日。
本文代码:https://github.com/zwdnet/2019-nCov-SIRmodel
与[7]作者讨论,我的算法是将S0与β作为独立的两个变量用两个循环嵌套分别遍历,他的做法是用每个S0的值代入微分方程算出相应的β值。他的算法应该更好一些,我正在尝试。另外在微信公众号上看到一篇更系统的关于此次疫情的数学模型的文章:https://mp.weixin.qq.com/s/rgaJtA4jioLOCHs_oCauDg
再次声明:本文只是我个人在家无聊的游戏作品,不是正儿八经的预测。我也不是流行病学专业人士。祝疫情早日结束!武汉加油!中国加油!
总结
以上所述是小编给大家介绍的python模拟预测一下新型冠状病毒肺炎的数据,希望对大家有所帮助!
来源:https://www.cnblogs.com/zwdnet/p/12246646.html


猜你喜欢
- 本文实例讲述了Selenium元素的常用操作方法。分享给大家供大家参考,具体如下:Selenium是一个用于Web应用程序测试的工具。Sel
- 目录一个可能你似曾相识的场景提升办公效率的法宝编码实现谷歌翻译爬虫Python 操作系统剪贴板项目地址妈妈再也不用担心我的英语了。一个可能你
- 直接通过element-ui自带的上传组件结合js即可,代码如下:HTML: &l
- 导读我们在使用selenium打开google浏览器的时候,默认打开的是一个新的浏览器窗口,而且里面不带有任何的浏览器缓存信息。当我们想要爬
- 我一直建议每个开发者都要有写博客记笔记的习惯,一来可以沉淀知识,二来可以帮助别人,我使用过很多博客平台,也用Python开发过博客系统,就这
- 使用 Beanstalkd 作为消息队列服务,然后结合 Python 的装饰器语法实现一个简单的异步任务处理工具.最终效果定义任务:from
- 本文实例讲述了PHP常用函数之获取汉字首字母功能。分享给大家供大家参考,具体如下://获取汉字的首字母function getFirstCh
- 1.在本地建立项目 可使用Eclipse,Idea等开发工具创建项目 打开根目录到所在在工程名的下一级2.使用git 客户端 进
- 本文介绍了tensorflow中next_batch的具体使用,分享给大家,具体如下:此处给出了几种不同的next_batch方法,该文章只
- 原理:使用js的定时任务函数setInterval设置时间,然后触发关闭事件参数说明title:提示框的标题msg:提示信息内容ico:显示
- 文章介绍了flask框架中的cookie和session。Session是在服务器端保存的一个数据结构,用来跟踪用户的状态,这个数据可以保存
- 写完调用天气接口的demo之后,小程序调用天气接口并且渲染在页面,顺便再调用了一下美图的接口API:美图APIurlwxml:<vie
- 今天给大家分享的主题是用百度的接口实现图片的文字识别。1.环境和配置要求整体是用Python实现,所需要使用的第三方库包括aip、PIL、k
- 随机生成10位数密码,字母和数字组合import string>>> import random>>>
- 前言嗨,彦祖们,不会过圣诞了还是一个人吧?今天我们来讲一下如何用python来画一个圣诞树,学会就快给那个她发过去吧,我的朋友圈已经让圣诞树
- 1.实现效果2.实现原理echarts官网:series-lines注意:流动特效只支持非平滑曲线(smooth:false)series-
- 实际应用中,会遇到需要把表的某些行转换成列,或者把列转换成行的情况。比如一张表在数据库中是这样的:图1但是,需要的结果可能是这样:图2这个时
- Doug Bowman,Google的Visual Design Lead离职了,一封带有感 * 彩的离职信惹发了大家不少的讨论。甚至还有人用
- 想必大家都知道MSSQL中SA权限是什么,可以说是至高无上。今天我就它的危害再谈点儿,我所讲的是配合NBSI上传功能得到WebShell。在
- 代码很简洁,也很简单,就不多废话了。/** * 去除多余的0 */ function del0($s)&