疫情模型

前几日做的是关于第二天新增确诊病例的预测,但是发现那个实在是太简单了,因为疑似病例中最后被确诊成肺炎的比例基本是不变的,在excel上把数据点画上去,然后画一条最优拟合线就完事了。

昨天突然意识到,自己其实有能力做一个关于疫情整体情况的模型,所以今天就动手操作一下。

模型

这次模型中有五类人:

    S- suseptible(可感人群)

    I- infectious(可以感染他人,但是由于种种原因没有被隔离,比如处于潜伏期,或者病情较轻,没有放在心上等。没有任何官方的数据来表明这类人的数量)

    D - diagnosed(被确诊而且被隔离,无法传染其他人。官方对于这类人的数量有明确的数据)

    R - recover(病情好转出院,获得免疫能力,无法被再次感染)

    DT - death(由于病重,不幸死亡)

有四个变化率:

    \alpha I -  每天被感染的人数

    \beta I - 每天被隔离治疗的人数

    \gamma D - 每天病情好转出院的人数

    \epsilon D - 每天死亡人数

接下来,我们可以写出他们变化率的表达式:

这一次我们只关注ID这两个量,不去计算疫情结束时的死亡和康复人数(如果想知道这两个量的话,看看死亡率和康复率的相对大小就知道了 ------ 人最后要么康复,要么死亡。按照现在的数据,如果不及时开发出新的治疗方法的话,后果真的很严重)

参数

好,接下来我们来看看如何确定这四个参数。

最好估计的就是\gamma (康复率)和\varepsilon (死亡率)。根据丁香园的数据,

\varepsilon  = 106/4565 = 0.0232

\gamma =62/4565 = 0.0135

对于\alpha\beta,我们可以做以下推导:

这里我们可以得到的量只有D_{t+1}D_{t}。也就是说,对于相同的CD_{t-1}-D_{t}\beta ,\alphaI_{t-1}可以有不同的解。不同的参数对应的模型肯定不一样,但是好在不是所有的模型都是最好的,我们可以想办法通过先找到最符合现实的模型,反过来求出这些参数值,然后使用通过这种方法求出的参数来进行预测。

思路:

\alpha,\beta和起始的I必须先确定其中的两个,才能求第三个值。那么先确定哪个呢?

\alpha和\beta是最难估计的了,那么我们就先确定初始的感染人数吧。新型冠状病毒的潜伏期的平均值有的专家说是10天,有的说是7天,我们选择保守一点的估计,7天。也就是说:

                                                                I_{initial} = D_{initial  }\div \frac{1}{7}

第二个参数我们选择\beta\beta的取值范围应该在0-1之间,如果大于1,那么就是说同一批人在感染了另一批人之后的同一天住进了医院,这明显不符合事实,那么我们就可以用暴力的方法,将0-1的值遍历一遍,每一遍都求一个最优的\alpha,最后看看当\beta取哪个值的时候这三个参数是最优的

公式推导:

D关于\alpha的导数:


设置损失函数并对\alpha求导:

\alpha更新公式:

其中theta是学习率(learning rate)  

代码实现:


模拟每天的变化:

计算损失:

更新\alpha

数据:

在建模的过程中最令我头疼的的就是数据了,搜遍了全网也也没有发现任何一家媒体把疫情的历史数据整理公开,虽然丁香园提供了历史数据图,但是并没有把所有的数据公开......

只有国家卫生健康委员会和央视新闻从1月19号开始每天播报确诊人数(可以从丁香园的实时播报查到:https://3g.dxy.cn/newh5/view/pneumonia_timeline?whichFrom=dxy)。也就是说,1月19号之前的数据要么根本就没有,要么就不可信......疫情可是从去年12月就开始了呀。

不管了,姑且先用这些数据吧


测试


初始化:

第一天(1月19号)的D(确诊人数)是291,假设潜伏期为7天,那么I就是2037。I和 \beta的取值先随便取,后面还要遍历搜索

遍历搜索:

看来当初始I = 2097,D = 291, 且\beta=0.02的时候,模型输出和真实数据拟合的是最好的。接下来把最好的\alpha也求出来吧。(在运行代码之前还要再初始化一遍)


我们知道,当\alpha - \beta大于0的时候,确诊人数会成指数上涨,现在看来情况依然不是很好啊

把最后拟合出来的图画出来康康:


嗯,看起来拟合的不错,模型的可信度还是有一些的

由于这世界上还没有没有任何关于I的数据,让我们用建好的模型来估计一下I的大小和历史变化吧


注意:这个结果是基于:模型结构没有问题,数学推导没有问题,数据没有问题,各种假设没有问题,初始人数没有问题等等才得到的。因此这张图仅供参考,本人无意制造恐慌。

根据模型输出的结果,截止到1月26号,已经有5万多人被感染了,只是这5万多人由于各种各样的原因没有被检测出来。大家要注意保护好自己呀,现在被确诊的人数可能仅仅只是所有被感染人群的冰山一角

至于什么时候情况会好转呢?现在从模型输出的数据上来看还没有好转的迹象,但是从丁香园的数据来看,疑似人数的增长已经放缓了

而且钟老已经放出话了:

也许明天就是转折点呢!


--------------------------------------------------2月5日更新

上面的求导和代码有些错误:

1)求导的时候少乘了个负1

2)代码也写的有问题,改正过代码在这里:


最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容

  • 前几日做的是关于第二天新增确诊病例的预测,但是发现那个实在是太简单了,因为疑似病例中最后被确诊成肺炎的比例基本是不...
    Rich070阅读 1,573评论 0 0
  • 首页 资讯 文章 资源 小组 相亲 登录 注册 首页 最新文章 IT 职场 前端 后端 移动端 数据库 运维 其他...
    Helen_Cat阅读 9,342评论 1 10
  • 点击查看原文 Web SDK 开发手册 SDK 概述 网易云信 SDK 为 Web 应用提供一个完善的 IM 系统...
    layjoy阅读 14,731评论 0 15
  • 很孤单很孤单很孤单。想说话又不想说话,不想说话又想说话。好像有你又好像没有你,没有你又好像有你。没有你孤单,有你更...
    大兔鱼儿阅读 2,444评论 1 2
  • 【导读】眼睛跟其他部位一样,也会时不时出点小毛病。一旦涉及到眼睛,所有问题都成了大问题。为了帮助每个人安全用药...
    爱眼小可爱阅读 5,792评论 0 11