前言

湖北新增一万多例中,大多数是临床病例,而真正的确诊病例,却只有两千多一点,比昨天略低。(关于临床病例和确诊病例的区别,请点击这里。)

不过不论如何,确诊人数数据由于临时加上临床诊断病例,导致目前没法用来做模型预测。现在能用来做预测的,只有 疑似病例 人数数据。

下面,我们将用疑似病例数据来做预测。原因如下:

  • 临床诊断将计入确认人数,这就意味着确认病例人数受到干扰,而且干扰的程度非常大~
  • 确诊病例来源于疑似病例。在《新型冠状病毒感染的肺炎诊疗方案(试行第五版)》中,在湖北省的病例诊断分类中增加了“临床诊断”,临床诊断来源于疑似病例,确诊病例来源于临床诊断和疑似病例。
  • 确诊病例在诊断过程受到硬件条件(如核酸试剂的数量、其他医疗设备的数量)的限制,遗漏的人数应该较多;而疑似病例判断标准则较低,只要其流行病史和是否发热等,遗漏的人数应该较少。

因此,我们将会用疑似病例来进行后续预测。

步骤一:抓取数据

从腾讯实时疫情数据,我们抓取 每日增加 数据,抓取结果如下:

    date  confirm  suspect  dead  heal  deadRate  healRate
0    1.2       77       27     0     0       0.0       0.0
1   1.21      149       53     3     0       2.0       0.0
2   1.22      131      257     8     0       6.1       0.0
3   1.23      259      680     8     6       3.1       2.3
4   1.24      444     1118    16     3       3.6       0.7
5   1.25      688     1309    15    11       2.2       1.6
6   1.26      769     3806    24     2       3.1       0.3
7   1.27     1771     2077    26     9       1.5       0.5
8   1.28     1459     3248    26    43       1.8       2.9
9   1.29     1737     4148    38    21       2.2       1.2
10   1.3     1982     4812    43    47       2.2       2.4
11  1.31     2102     5019    46    72       2.2       3.4
12  2.01     2590     4562    45    85       1.7       3.3
13  2.02     2829     5173    57   147       2.0       5.2
14  2.03     3235     5072    64   157       2.0       4.9
15  2.04     3893     3971    65   262       1.7       6.7
16  2.05     3697     5328    73   261       2.0       7.1
17  2.06     3143     4833    73   387       2.3      12.3
18  2.07     3401     4214    86   510       2.5      15.0
19  2.08     2656     3916    89   600       3.4      22.6
20  2.09     3062     4008    97   632       3.2      20.6
21   2.1     2484     3536   108   716       4.3      28.8
22  2.11     2022     3342    97   744       4.8      36.8

附抓取代码(类似过程可见《新型冠状病毒数据抓取及整理详细流程》):

import json
import pandas as pd
import numpy as np
from datetime import datetime
# 这里requests-html需要python3.5以上版本
from requests_html import HTMLSession# 创建会话
session = HTMLSession();
url = r'https://view.inews.qq.com/g2/getOnsInfo?name=disease_h5';
# 从目标网址中读取数据
r = session.get(url).json();
# 将数据存储为json格式
data = json.loads(r['data']);
# 每日新增数据(将'date'列排在第一列)
chinaDayAddData = pd.DataFrame(data['chinaDayAddList'])[['date', 'confirm', 'suspect', 'dead', 'heal', 'deadRate', 'healRate']];
# 保存
chinaDayAddData.to_csv(r'F:\开课吧\RS基础课\2019_CoV_chinaDayAddData.csv', encoding='gbk', index=None)
# 读取
chinaDayAddData = pd.read_csv(r'F:\开课吧\RS基础课\2019_CoV_chinaDayAddData.csv', encoding='gbk');
chinaDayAddData['date'] = chinaDayAddData['date'].map(str);

步骤二:估计实际确诊人数

我们要做以下事情:

  • 估计疑似病例中确诊的比例
  • 根据每日增加数据确定确诊人数
  • 根据SI模型,确定最终确诊人数、疫情结束的大概时间

我们首先看如何通过现有数据,估计疑似病例中确诊比例。

# 最小二乘法
# 采样点
# 确诊病例每日增加数
Yi = np.array(list(chinaDayAddData['confirm']));
# 疑似病例每日增加数
Xi = np.array(list(chinaDayAddData['suspect']));# 需要拟合的函数func及误差error
def func(p, x):return p * x;def error(p, x, y):return func(p, x) - y;# 疑似病例确诊比例的初始值
p0 = 0.5;
# 最小二乘法做出的疑似转确诊比例
percent = leastsq(error, p0, args=(Xi, Yi))[0][0];
print(percent);

这样得到的疑似转确诊比例为0.5960971132108407。

其次,我们根据每日增加的每日疑似病例增加数据确定确诊人数。

  • 确定每日疑似病例总数,我们需要累加和函数cumsum()
  • 确定每日确诊病例总数,我们只需将0.5960971132108407乘以每日疑似病例总数

代码如下:

predict = chinaDayAddData['suspect'].cumsum() * 0.6

我们将现在作出的预测与官方给出的每日确诊人数进行比较,如下:

可见,我们预测的每日确诊病例总数要比官方给出的要高,但是最近几天开始重合。这其实可以这样理解,在当前国家高度重视和充分宣传之下,个人和各级政府不再存在瞒报,因此可以认为是真实数据。

退一万步,如果所有疑似病例都转变成确诊病例,那么如绿线所示,当前的实际感染总人数为7万多人~

步骤三:预测

这里,我们选择SI模型。在说明选择这个最简单模型之前,先给一个简介。

SI模型

在这个最简单的传染病模型中,我们仅仅考虑健康人(即易感者,Susceptibles)和病人(即感染者,Infection)。

我们作出如下假设:

  • 总人数NNN不变,病人和健康人的比例分别为i(t)i(t)i(t)和s(t)s(t)s(t)。其中,ttt为从0开始的时间,以天为单位。
  • 每个病人每天接触到的人数为λ\lambdaλ,且能使得其中的健康人感染。

根据上述假设,我们可以列出下式:
N⋅[i(t+Δt)−i(t)]=λs(t)⋅Ni(t)ΔtN\cdot [i(t+\Delta t)-i(t)]=\lambda s(t) \cdot Ni(t)\Delta tN⋅[i(t+Δt)−i(t)]=λs(t)⋅Ni(t)Δt
经过简单的变形,我们有
didt=λsi\frac{d i}{d t}=\lambda s idtdi​=λsi
且根据假设,
i+s=1i+s=1i+s=1
因此,我们可以写出微分方程,如下
{didt=λi(1−i)i(0)=i0\left\{ \begin{array}{lll} \frac{d i}{d t}&=&\lambda i(1-i)\\ i(0)&=&i_0 \end{array}\right. {dtdi​i(0)​==​λi(1−i)i0​​
其中,i(0)=i0i(0)=i_0i(0)=i0​为初始条件。

求解上述微分方程,可以得到Logistic模型,即
i(t)=11+(1i0−1)e−λt=i0eλti0(eλt−1)+1\begin{array}{lll} i(t)&=&\frac{1}{1+\left(\frac{1}{i_0}-1\right)e^{-\lambda t}}\\ &=&\frac{i_0e^{\lambda t}}{i_0(e^{\lambda t}-1)+1} \end{array} i(t)​==​1+(i0​1​−1)e−λt1​i0​(eλt−1)+1i0​eλt​​

但通过观察这个公式,我们发现,随着时间ttt的增加,感染比例i(t)i(t)i(t)也逐渐增加,最后趋近于1。这就意味着最后所有人都会被感染。

说明

可以反向思考,在疫情结束之后看现在,有什么特点呢?

  • 假设整个疫情造成的确诊人数为NNN;
  • 当前这NNN个人中,我们假设确诊比例为i(t)i(t)i(t),尚且健康的人的比例为s(t)s(t)s(t);
  • 在这一过程中,确诊比例不断增大,且存在斜率先增加后减少的过程;
  • 最终i=1,s=0i=1, s=0i=1,s=0。

SI模型预测

我们用SI模型根据估计的确诊人数进行预测。我们需要估计三个参数N,i0,λN, i_0, \lambdaN,i0​,λ。

# 采样点
Yi = np.array(predict);
Xi = np.arange(len(Yi));# Logistic模型
def i_t(t, i0, lamb):i1 = i0 * math.exp(lamb * t);return i1 / (1 + i1 + i0);# 需要拟合的函数func及误差error
def func(p, x):i0, N, lamb = p;y = np.array([i_t(t, i0, lamb) * N for t in x]);return ydef error(p, x, y):return func(p, x) - y;p0 = [0.01, 1000, 1];
#print(func(p0, Xi))
#print(error(p0, Xi, Yi));
Para = leastsq(error, p0, args=(Xi,Yi));
print(Para);

这样得到的值为(i0,N,λ)=(1.88087189e−02,4.72770019e+04,2.90671789e−01)(i_0, N, \lambda)=(1.88087189e-02, 4.72770019e+04, 2.90671789e-01)(i0​,N,λ)=(1.88087189e−02,4.72770019e+04,2.90671789e−01)。

将这样的值代入原方程,我们可以作出

# 预测1月20号之后40天的确诊人数的状况
Xi = np.arange(40);
y_fit = func(Para[0], Xi);
plt.plot(Xi, y_fit, label=r'predict_confirm');
#plt.plot(Xi, Yi, label=r'original');
plt.legend();
plt.show();


即,在3月中旬左右疫情不会再大量增加,同时,最终确诊总人数为47277人。

同样,我们退一万步说,疑似病例全部确诊为确诊病例,这样,我们有


在极端情况(所有疑似病例最终都确诊为新冠肺炎)下,最终确诊总人数为78744,同样在3月中旬左右结束。

新冠肺炎疫情数学模型的一点想法相关推荐

  1. 居家学习:新冠肺炎疫情下中国高校基于直播的远程教育体验的混合方法分析

    Learning from Home: A Mixed-Methods Analysis of Live Streaming Based Remote Education Experience in ...

  2. 新冠肺炎疫情把科研推上“云端”

    来源:新华网 美国威斯康星国家灵长类动物研究中心的戴夫·奥康纳清晨收到在伦敦的一名合作伙伴发来的论文预印本.这项研究在中国完成,两人通过企业协同云端办公软件Slack讨论了一上午.下午2点,奥康纳打开 ...

  3. 新冠肺炎疫情预测与防控策略评价

    新冠肺炎疫情预测与防控策略评价 新冠肺炎疫情在全球蔓延以来,各个国家采取了不同的防疫策略和措施, 导致新冠肺炎疫情感染人数各异,效果不尽相同. 问题1:请你们团队收集几个典型国家(如:美国.英国.印度 ...

  4. 含最新数据! 使用Python检测新冠肺炎疫情拐点

    注:本文案例仅供技术学习,不代表研究性观点. 本文对应代码.数据及文献资料已上传至Github仓库https://github.com/CNFeffery/DataScienceStudyNotes ...

  5. 每日一练:Python爬虫爬取全国新冠肺炎疫情数据实例详解,使用beautifulsoup4库实现

    Python 爬虫篇 - 爬取全国新冠肺炎疫情数据实例详解 效果图展示 第一章:疫情信息的下载与数据提取 ① 爬取页面数据到本地 ② json 字符串正则表达式分析 ③ 提取数据中的 json 字符串 ...

  6. 第四范式荣获“工业和信息化系统抗击新冠肺炎疫情先进集体”称号

    近日,在全国工业和信息化工作会议上,工信部党组为工业和信息化系统抗击新冠肺炎疫情工作中,对疫情防控.促进经济持续稳定恢复发挥重要作用的全国100家「先进集体」和300名「先进个人」进行了表彰,凭借在疫 ...

  7. 最新!兰州大学发布对上海市的新冠肺炎疫情预测!

    这段时间,上海市疫情牵动着所有人的心.据数据显示,自 2022 年 3 月 1 日上海市报告新冠肺炎本土确诊病例和本土无症状感染者以来,截至 2022 年 4 月 10 日 24 时,上海市已累计报告 ...

  8. 视频分享|新冠肺炎疫情趋势预测大赛技术详解

    点击我爱计算机视觉标星,更快获取CVML新技术 本篇文章的内容是???? 3月1号FlyAI直播视频整理分享 视频内容: ????新冠肺炎疫情趋势预测大赛技术详解 ????用户问题答疑 前半部分会带领 ...

  9. 谷歌设立 COVID-19资助金,助力漏洞猎人度过新冠肺炎疫情

     聚焦源代码安全,网罗国内外最新资讯! 编译:奇安信代码卫士团队 本周,谷歌宣布设立 COVID-19 资助金,为参与其漏洞奖励计划 (VRP) 的漏洞猎人提供帮助. 这项新推出的资助金是对漏洞研究资 ...

最新文章

  1. 环形、扭环、LFSR计数器
  2. PHP 安装 扩展时 抛出 /usr/local/Cellar/php@7.1/7.1.25/pecl 异常解决
  3. HALCON示例程序gray_features.hdev提取灰度图的不同特征(area_center_gray 、elliptic_axis_gray)
  4. .net mysql字符串截取_【MySQL】字符串截取之SUBSTRING_INDEX和【MySQL】字符串四则运算...
  5. C#获取文件(磁盘驱动器)的关联图标(使用API SHGetFileInfo)
  6. 构建之法首周阅读体会
  7. 玩转SpringBoot 2 之项目启动篇
  8. 尼尔机械纪元机器人驱动_中国的工业机器人发展到了哪一步?
  9. leetcode 之Rotate List(18)
  10. 应用架构、业务架构、技术架构和业务流程图详解
  11. 从计算机网络系统组成看 计算机网络可分为,从计算机网络系统组成的角度看,计算机网络可以分为 子网和资源子网。...
  12. 什么软件测试情侣头像,情侣头像搜索另一半 另一半情头查找器在线匹配 - 株洲娱乐网...
  13. win10电脑双屏如何设置不同的桌面
  14. $body.on(‘click‘,‘button‘,function() {中.on在这里是什么意思的
  15. 单元测试报错Wanted but not invoked:
  16. 新的一年,知足,上进,不负野心!
  17. php 图片上载 wordpress,WordPress 使用 Jcorp上传并裁剪图片作为自定义头像 —— PHP 后台部分...
  18. 硬件使用74hc138的C语言程序,【Arduino教程】第三十一讲:74HC138实验
  19. 机械振动信号15个时域统计指标
  20. “向前兼容”与“向后兼容”

热门文章

  1. 华硕和超微服务器主板稳定性,华硕、超微、泰安三款双路服务器主板横测
  2. radius java_Java处理Radius access-challenge
  3. 度度熊的01世界 DFS
  4. Aizu ITP2_6_A 二分
  5. C++中1LL避免强制转换
  6. 【转】Java 8十个lambda表达式案例
  7. 【B/S】ASP.NET---验证码-生成和切换
  8. Log4J基础详解及示例大全
  9. 生活大爆炸系列之制作望远镜架
  10. python-基于UDP通信的套接字,socketserver模块的使用