用python实现传染病模型

1.SI模型
- 1.1 代码实现
- 1.2 模型的结果
2.SIS （治愈后仍然还是易感者）
- 2.1 代码实现
- 2.2模型的结果
3 SIR模型（治愈后直接移除）
- 3.2代码实现
- 3.2绘制图像：
4.SEIR 模型（新增一个人群，叫潜伏者E）
- 4.1代码实现
- 4.2模型的结果

参考文章为：https://www.kesci.com/mw/project/60161bf6ac79f40016b7d7d9

1.SI模型

示意图：

我们假设城市有一千万（N=10的7次方）人，每个患者每天接触感染每天0.8人（lamda=0.8），初始感染人数为45人(i0 = 45/N)，我们来模拟70天(T=70)的情况。

1.1 代码实现

import numpy as np
import matplotlib.pyplot as plt
#population N = 1e7#simuation TimeT = 70
# susceptiable ratio 易感染比例
s = np.zeros([T])
# infective ratio 感染比例
i = np.zeros([T])# 由我们的实验可以知道，lamda是有取值范围的，似乎超过1就不可以了
# concat ratelamda = 0.9# initial infective people i[0] = 45.0 / N
s[0] = 1-i[0]for t in range(T-1):i[t+1] = i[t] + i[t]  * lamda * (1.0 - i[t])s[t+1] = 1 - i[t+1]

1.2 模型的结果

#感染者随着天数的变化曲线
fig, ax = plt.subplots(figsize=(8,4))
ax.plot(i, c='r', lw=2)
ax.plot(s, c='b', lw=2)ax.set_xlabel('Day',fontsize=20)
ax.set_ylabel('Infective Ratio', fontsize=20)
ax.grid(1)
plt.xticks(fontsize=20)plt.yticks(fontsize=20);

由上图可见，大约在25天左右，全部人群都会变成感染者，感染率为1

lamda的现实意义就是该城市的卫生水平，衡量的是消毒，隔离这些措施执行得怎么样。

回到传染病模型，按照SI模型计算的结果，我们全人类都会患病，这好可怕！原因是我们忽略了一个很重要的因素，那就是我们有奋斗在一线的医护人员，我们会被治愈！所以SI模型只适合研究具有高传染风险又不能被治愈的病（比如HIV）。

但是对于其他病，我们是可以靠医疗和自身免疫系统康复的，那么紧接着的一个问题就是，被治愈后还会再被传染上嘛？根据这个问题的回答不同，我们有了两个不同的模型，SIR 和 SIS。现在可以揭晓，SIR的R的含义了，就是移出者（Removed），现实含义就是指被治愈后不会再被感染的人。而SIS表示治愈后仍然还是易感者。下面我们用python来分别实现这两个模型。

2.SIS （治愈后仍然还是易感者）

为了实现这个模型，我们需要引入新的一个参数，治愈率γ\gammaγ。好啦，先上我们的新示意图：

和SI模型做比较，区别就是计算感染者的增加数时要减去被治愈的人数。
所以这时候每天的增加的感染者为：λ×iN×s−γ×iN\lambda \times i N \times s-\gamma \times i Nλ×iN×s−γ×iN ,
增加的感染率为：λ×i×s−γi\lambda \times i \times s-\gamma iλ×i×s−γi 。
模型完成啦，修改python代码：

2.1 代码实现

# susceptiable ratio
s = np.zeros([T])# infective ratioi = np.zeros([T])# concat rate(这个其实是有效的感染率)lamda = 1.0# recover rate(治愈率)gamma = 0.5# initial infective people i[0] = 45.0 / N
s[0] = 1-i[0]for t in range(T-1):i[t+1] = i[t] + i[t] * (1- i[t])* lamda - i[t] * gammas[t+1] = 1 - i[t+1]

2.2模型的结果

 fig, ax = plt.subplots(figsize=(8,4))
ax.plot(i, c='r', lw=2)
ax.plot(s, c='b', lw=2)ax.set_xlabel('Day',fontsize=20)
ax.set_ylabel('Infective Ratio', fontsize=20)
ax.grid(1)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20);

可以看到，达到最大感染率的时间退后10天左右，最后感染和治愈达到动态平衡，人群中有始终有一半的人感染着。所以，SIS模型适合研究具有传染性和反复性的流行病，比如常见流感。同样的，感兴趣的话，改变lamda和gamma的值，观察曲线的变化。和lamda不同的是，gamma的现实意义就是对这种疾病的治疗水平。

3 SIR模型（治愈后直接移除）

加入了移出者，被治愈的病人不会再被传染，先上我们的新示意图：

SIR 模型
注意到这里，人群被分成了三类，不再只有I和S，所以相比于之前的模型，我们需要找到新的约束关系。现在我们需要分别计算三种人每天的增加量了：

易感者：每天都在被传染，所以一直在减少，减少量为被传染的人数：λNis\lambda N i sλNis
感染者：增加了被感染的人，减少了治愈的人： λNis−γNi\lambda N i s-\gamma N iλNis−γNi
移出者：增加了治愈的人： γNi\gamma N iγNi

建模完成，修改python代码，并且假设人群普遍易感，新型疾病，初始没有移出者。

γ=0.0821,λ=0.2586，初始易感人数为一千万,初始感染10人，那么我们的城市总人数N=1e7+10\gamma=0.0821,\lambda=0.2586 ，初始易感人数为一千万, 初始感染10人，那么我们的城市总人数 N=1 e 7+10γ=0.0821,λ=0.2586，初始易感人数为一千万,初始感染10人，那么我们的城市总人数N=1e7+10

3.2代码实现

# populationN = 1e7 + 10 # simuation Time / Day
T = 170# susceptiable ratios = np.zeros([T])# infective ratioi = np.zeros([T])#remove ratio
r = np.zeros([T])# contact ratelamda = 0.2586# recover rate gamma = 0.0821# initial infective peoplei[0] = 10.0 / Ns[0] = 1e7 / Nfor t in range(T-1):i[t+1] = i[t]+i[t]* lamda * s[t] - gamma*i[t]s[t+1] = s[t] - lamda * i[t]*s[t]r[t+1] =r[t] + gamma * i[t]

3.2绘制图像：

fig, ax = plt.subplots(figsize=(10,6))
ax.plot(s, c='b', lw=2, label='S')
ax.plot(i, c='r', lw=2, label='I')
ax.plot(r, c='g', lw=2, label='R')
ax.set_xlabel('Day',fontsize=20)
ax.set_ylabel('Infective Ratio', fontsize=20)
ax.grid(1)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.legend();

感染人数峰值发生在一个月左右，最大感染人数不到人群的20%，但是最终人群的80%都会得此病（就是最终的移出者的比例）。SIR模型适合研究没有潜伏期的急性传染病，治疗后能够痊愈并具有抗病性。

4.SEIR 模型（新增一个人群，叫潜伏者E）

但是，SIR模型和实际情况的出入会比较大，因为忽略了太多因素了，比如说潜伏期，比如说政策调控，药物，出生死亡等等。下面我们可以和前面一样，把潜伏期考虑进去，新增一个人群，叫潜伏者E（exposed）：

SEIR模型
同样的我们需要计算各人群每天的增加量：

S：每天减少： λNis\lambda N i sλNis

E：每天增加传染，减少发病： λNis−σNe\lambda N i s-\sigma N eλNis−σNe

I：每天增加发病，减少治愈： σNe−γNi\sigma N e-\gamma N iσNe−γNi

R：每天增加治愈： γNi\gamma N iγNi

建模完成，修改我们的python程序，这里的 σ\sigmaσ可以理解为潜伏期的倒数。给的4天。新型冠状病毒给目前临床的潜伏期是3-14天。

4.1代码实现

# population
N = 1e7 + 10 + 5
# simuation Time / Day
T = 170
# susceptiable ratio
s = np.zeros([T])
# exposed ratio
e = np.zeros([T])
# infective ratio
i = np.zeros([T])
# remove ratio
r = np.zeros([T])# contact rate
lamda = 0.5
# recover rate
gamma = 0.0821# 潜伏期# exposed periodsigma = 1 / 4# initial infective people
i[0] = 10.0 / N
s[0] = 1e7 / N
e[0] = 40.0 / Nfor t in range(T-1):s[t + 1] = s[t] - lamda * s[t] * i[t]e[t + 1] = e[t] + lamda * s[t] * i[t] - sigma * e[t]i[t + 1] = i[t] + sigma * e[t] - gamma * i[t]r[t + 1] = r[t] + gamma * i[t]

4.2模型的结果

fig, ax = plt.subplots(figsize=(10,6))
ax.plot(s, c='b', lw=2, label='S')
ax.plot(e, c='orange', lw=2, label='E')
ax.plot(i, c='r', lw=2, label='I')
ax.plot(r, c='g', lw=2, label='R')
ax.set_xlabel('Day',fontsize=20)
ax.set_ylabel('Infective Ratio', fontsize=20)
ax.grid(1)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.legend();

用python实现传染病模型传染病模型相关推荐

vecm模型怎么写系数_经典传染病的SIR模型（基于MATLAB）
经典的SIR模型是一种发明于上个世纪早期的经典传染病模型,此模型能够较为粗略地展示出一种传染病的发病到结束的过程,其核心在于微分方程,本次我们利用matlab来对此方程进行其中三个主要量 S是易感人 ...
Python基于statsmodels包构建多元线性回归模型：模型构建、模型解析、模型推理预测
Python基于statsmodels包构建多元线性回归模型:模型构建.模型解析.模型推理预测目录
python可视化多个机器学习模型在独立测试集（test data set）上面的AUC值、可视化模型效能
python可视化多个机器学习模型在独立测试集(test data set)上面的AUC值.可视化模型效能 # x_lables为模型名称列表,包括,逻辑回归.朴素贝叶斯.支持向量机.随机森林.xgb ...
python可视化多个机器学习模型在训练集(train set)上交叉验证（cross validation）的AUC值、可视化模型效能
python可视化多个机器学习模型在训练集(train set)上交叉验证(cross validation)的AUC值.可视化模型效能 # 所有的模型中填写的参数都是通过randomsearchcv ...
python学习笔记（IO模型）
1.IO模型介绍: io模型一般有五种: * blocking IO * nonblocking IO * IO multiplexing * s ...
python多线程实现生产者消费者_用Python实现多线程“生产者-消费者”模型的简单例子...
用 Python 实现多线程"生产者 - 消费者"模型的简单例子生产者消费者问题是一个著名的线程同步问题, 该问题描述如下: 有一个生产者在生产产品, 这些产品将提供给若干个消费 ...
手把手教程：用Python开发一个自然语言处理模型，并用Flask进行部署
截住到目前为止,我们已经开发了许多机器学习模型,对测试数据进行了数值预测,并测试了结果.实际上,生成预测只是机器学习项目的一部分,尽管它是我认为最重要的部分.今天我们来创建一个用于文档分类.垃圾过滤的 ...
python回归分析预测模型_在Python中如何使用Keras模型对分类、回归进行预测
姓名:代良全学号:13020199007 转载自:https://www.jianshu.com/p/83ba11abdffc [嵌牛导读]: 在Python中如何使用Keras模型对分类.回归进行 ...
python加载模型_解决python 无法加载downsample模型的问题
downsample 在最新版本里面修改了位置 from theano.tensor.single import downsample (旧版本) 上面以上的的import会有error raise: ...
python信用评分卡_基于Python的信用评分卡模型分析（二）
上一篇文章基于Python的信用评分卡模型分析(一)已经介绍了信用评分卡模型的数据预处理.探索性数据分析.变量分箱和变量选择等.接下来我们将继续讨论信用评分卡的模型实现和分析,信用评分的方法和自动评分 ...

用python实现传染病模型传染病模型