传染病发展趋势与疾控决策数值模拟

于 2024-12-18 10:30:00 发布
阅读量1.1k 收藏 47
点赞数 51
文章标签： python 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/m0_62222486/article/details/144491887
版权

第1关：传染病模型 SIR 模型-前向欧拉法求解

任务描述
新型冠状病毒感染的肺炎疫情引发全球广泛关注。2020 年 3 月 11 日，世界卫生组织总干事谭德塞宣布，根据评估，世卫组织认为当前新冠肺炎疫情可被称为全球大流行（ pandemic ）。我们可以利用传染病模型，对疾病的发展进行简单预测，这对疫情发展趋势分析具有一定的参考价值。

本关的目标就是让学习者利用 Python 编程实现经典传染病模型—— SIR模型（ Susceptible Infected Recovered Model ）对疫情趋势进行分析。本关中，我们采用前向欧拉法完善 SIR 模型，对疫情发展进行模拟。

编程要求
利用前向欧拉法完善 SIR 模型，实现 SIR 类，其属性

beta：对应SIR模型中的 β；
gamma：对应SIR模型中的γ；
y0：对应SIR模型中的(S, I, R)人数。
并编写方法：

f(self, y): 实现 SIR 模型，其中参数 y 是当前状态 [S, I, R] ，该函数需返回 [S’,I’,R’]；

solve(self, x): 实现前向欧拉法，其中参数 x 是列表 [0, …, t-1] ，该函数返回 t 天内每天感染人数的列表 (列表大小为天数t)。（注意：本场景中应用前向欧拉法时， forwardEuler() 函数第 5 行 f(x[k],y[k])`` 中 x[k]对f` 函数的计算没有影响）。

任务提示
测试输入：

N = 1e8 # 武汉总人数：1000万人
gamma = 1/25 # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0 = [N-1, 1, 0] # 初始发病1人，其他人员正常 [S0, I0, R0]
t = range(0,5,1) # 模拟5天的发展情况，单位时间为1天
beta = 1.0/N # 平均传染率
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y = simulation.solve(t)
预期输出：

[1.0000, 1.9600, 3.8416, 7.5295, 14.7579]

# import ode_7 as ode
import numpy as np

def myprint_list(l):
    if type(l)!=list:
        print('Error: the result is not a list')
    else:
        print('[', end='')
        for i in range(len(y_f)-1):
            print("%.4f, "%y_f[i], end='')
        print("%.4f]"%y_f[len(y_f)-1])
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self,y): # y为当前状态，即列表[S,I,R]
        # ------------------begin-----------------------
        y1=-self.beta*y[0]*y[1]
        y2=self.beta*y[0]*y[1]-self.gamma*y[1]
        y3=self.gamma*y[1]
        return np.array([y1,y2,y3])

        
        # -------------------end------------------------
    
    def solve(self, t):
        # ------------------begin-----------------------
        y = [self.y0] * len(t)
        for k in range(len(t) - 1):
            h = t[k + 1] - t[k]
            y[k + 1] = y[k] + h *self.f(y[k])
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res
        

        # -------------------end------------------------


N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗） 
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常, 即[S0, I0, R0]
beta = 1.0/N            # 感染系数
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
for T in [5, 10, 20]:
    t   = range(0, T, 1) # 模拟T天的发展情况，单位时间为1天
    y_f = simulation.solve(t)
    myprint_list(y_f)

第2关：传染病模型 SIR 模型-后向欧拉法求解

任务描述
新型冠状病毒感染的肺炎疫情引发全球广泛关注。2020 年 3 月 11 日，世界卫生组织总干事谭德塞宣布，根据评估，世卫组织认为当前新冠肺炎疫情可被称为全球大流行（ pandemic ）。我们可以利用传染病模型，对疾病的发展进行简单预测，这对疫情发展趋势分析具有一定的参考价值。

测试输入：

N = 1e8 # 武汉总人数：1000万人
gamma = 1/25 # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0 = (N-1, 1, 0) # 初始发病1人，其他人员正常 (S0, I0, R0)
t = range(0, 5, 1) # 模拟5天的发展情况，单位时间为1天
beta = 1.0/N # 平均传染率
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y = simulation.solve(t)

import numpy as np

def myprint_list(l):
    if type(l)!=list:
        print('Error: the result is not a list')
    else:
        print('[', end='')
        for i in range(len(y_f)-1):
            print("%.4f, "%y_f[i], end='')
        print("%.4f]"%y_f[len(y_f)-1])
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self, y):  # y为当前状态，列表[S,I,R]
        # ------------------begin-----------------------
        y1=-self.beta*y[0]*y[1]
        y2=self.beta*y[0]*y[1]-self.gamma*y[1]
        y3=self.gamma*y[1]
        return np.array([y1,y2,y3])



        # -------------------end------------------------
    
    def solve(self, x):
        # ------------------begin-----------------------
        y = [self.y0] * len(t)
        for k in range(len(t)-1):
            h = t[k+1] - t[k]
            yp=y[k]+h*self.f(y[k])
            y[k+1] = y[k] + h *self.f(yp)
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res
        






        # -------------------end------------------------


N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗） 
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常, 即[S0, I0, R0]
beta = 1.0/N            # 平均传染率倒数
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
for T in [5, 10, 20]:
    t   = range(0, T, 1) # 模拟T天的发展情况，单位时间为1天
    y_f = simulation.solve(t)
    myprint_list(y_f)

第3关：传染病模型 SIR 模型-改进欧拉法(梯形法)求解

任务描述
新型冠状病毒感染的肺炎疫情引发全球广泛关注。2020 年 3 月 11 日，世界卫生组织总干事谭德塞宣布，根据评估，世卫组织认为当前新冠肺炎疫情可被称为全球大流行（ pandemic ）。我们可以利用传染病模型，对疾病的发展进行简单预测，这对疫情发展趋势分析具有一定的参考价值。

本关的目标就是让学习者利用 Python 编程实现经典传染病模型—— SIR 模型（ Susceptible Infected Recovered Model ）对疫情趋势进行分析。本关中，我们采用改进的欧拉法（梯形法，两步）完善 SIR 模型，对疫情发展进行模拟。

import numpy as np

def myprint_list(l):
    if type(l)!=list:
        print('Error: the result is not a list')
    else:
        print('[', end='')
        for i in range(len(y_f)-1):
            print("%.4f, "%y_f[i], end='')
        print("%.4f]"%y_f[len(y_f)-1])
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self, y): #  y为当前状态，列表[S,I,R]
        S = y[0]
        I = y[1]
        R = y[2]
        return np.array([-self.beta * S * I, self.beta * S * I - self.gamma * I, self.gamma * I])
    
    def solve(self, x):
        # ------------------begin-----------------------
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1] - x[k]
            yp=y[k]+h*self.f(y[k])
            yc=y[k]+h*self.f(yp)
            y[k+1] = 1/2*(yp+yc)
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res






        # -------------------end------------------------


N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗） 
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常 [S0, I0, R0]
t     = range(0, 11, 1) # 模拟10天的发展情况，单位时间为1天
beta = 1.0/N            # 平均传染率

simulation = SIR(beta=beta, gamma=gamma, y0=y0)
for T in [5, 10, 20]:
    t   = range(0, T, 1) # 模拟T天的发展情况，单位时间为1天
    y_f = simulation.solve(t)
    myprint_list(y_f)

第4关：引入疫情防控措施-隔离机制

任务描述
为有效防止疫情进一步扩散和蔓延，各地纷纷出台各种举措加强疫情防控工作。为了更准确的分析数据，预测疫情发展，本关将在第3关的基础上引入隔离机制，模拟疫情防控。

测试输入：

N = 1e8 # 武汉总人数：1000万人
gamma = 1/25 # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0 = [N-1, 1, 0] # 初始发病1人，其他人员正常 [S0, I0, R0]
beta = 1.0/N # 平均传染率
t = range(0, 5, 1) # 模拟5天的发展情况，单位时间为1天
print(“———————————— 无隔离下感染情况：————————————”)
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y_f = simulation.solve_with_quarantine(t)
myprint_list(y_f)
print(“———————————— 隔离确诊患者：————————————”)
gamma1 = 1/15 # 隔离确诊患者：按最长15天发病确诊后被隔离
y_f = simulation.solve_with_quarantine(t, gamma1)
myprint_list(y_f)
print(“———————————— 隔离疑似人员：————————————”)
gamma1 = 1/3 # 隔离疑似人员：按平均3天被隔离
y_f = simulation.solve_with_quarantine(t, gamma1)
myprint_list(y_f)
预期输出：

———————————— 无隔离下感染情况：————————————
[1.0000, 2.4208, 5.8603, 14.1865, 34.3428]
———————————— 隔离确诊患者：————————————
[1.0000, 2.3689, 5.6116, 13.2933, 31.4904]
———————————— 隔离疑似人员：————————————
[1.0000, 1.8889, 3.5679, 6.7394, 12.7299]

import numpy as np

def myprint_list(l):
    if type(l)!=list:
        print('Error: the result is not a list')
    else:
        print('[', end='')
        for i in range(len(y_f)-1):
            print("%.4f, "%y_f[i], end='')
        print("%.4f]"%y_f[len(y_f)-1])
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self, y): #  y为当前状态，列表[S,I,R]
        S = y[0]
        I = y[1]
        R = y[2]
        return np.array([-self.beta * S * I, self.beta * S * I - self.gamma * I, self.gamma * I])
    
    def solve_with_quarantine(self, x):
        # ------------------begin-----------------------
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1] - x[k]
            yp=y[k]+h*self.f(y[k])
            yc=y[k]+h*self.f(yp)
            y[k+1] = 1/2*(yp+yc)
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res





        # -------------------end------------------------

N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常 [S0, I0, R0]
beta = 1.0/N            # 平均传染率

for T in [5, 10, 20]:
    t   = range(0, T, 1) # 模拟T天的发展情况，单位时间为1天
    print("———————————— 无隔离下感染情况：————————————")
    simulation = SIR(beta=beta, gamma=gamma, y0=y0)
    y_f = simulation.solve_with_quarantine(t)
    myprint_list(y_f)
    print("———————————— 隔离确诊患者：————————————")
    gamma1 = 1/15            # 隔离确诊患者：按最长15天发病确诊后被隔离
    simulation = SIR(beta=beta, gamma=gamma1, y0=y0)
    y_f = simulation.solve_with_quarantine(t)
    myprint_list(y_f)
    print("———————————— 隔离疑似人员：————————————")
    gamma2 = 1/3            # 隔离疑似人员：按平均3天被隔离
    simulation = SIR(beta=beta, gamma=gamma2, y0=y0)
    y_f = simulation.solve_with_quarantine(t)
    myprint_list(y_f)

第5关：引入疫情防控措施-出行控制机制

任务描述
为有效防止疫情进一步扩散和蔓延，各地纷纷出台各种举措加强疫情防控工作。为了更准确的分析数据，预测疫情发展，本关将在第三关的基础上引入出行控制机制，模拟疫情防控。

测试输入：

N = 1e8 # 武汉总人数：1000万人
gamma = 1/25 # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0 = [N-1, 1, 0] # 初始发病1人，其他人员正常 [S0, I0, R0]
beta = 1.0/N # 平均传染率
t = range(0, 5, 1) # 模拟5天的发展情况，单位时间为1天
print(“—————————————— 无出行控制下感染情况：——————————————”)
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y_f = simulation.solve_with_control(t)
myprint_list(y_f)
for state in [2,5,10]: #出行控制的程度,对应传染率分别下降为原来的1/state
simulation = SIR(beta=beta, gamma=gamma, y0=y0)
y_f = simulation.solve_with_control(t,state)
print(“———出行控制后，传染率为原来的 1/{} 时,对应的感染情况：———”.format(state))
myprint_list(y_f)
预期输出：

—————————————— 无出行控制下感染情况：——————————————
[1.0000, 2.4208, 5.8603, 14.1865, 34.3428]
———出行控制后，传染率为原来的 1/2 时,对应的感染情况：———
[1.0000, 1.5658, 2.4517, 3.8389, 6.0110]
———出行控制后，传染率为原来的 1/5 时,对应的感染情况：———
[1.0000, 1.1728, 1.3755, 1.6131, 1.8919]
———出行控制后，传染率为原来的 1/10 时,对应的感染情况：———
[1.0000, 1.0618, 1.1274, 1.1971, 1.2711]

import numpy as np

def myprint_list(l):
    if type(l)!=list:
        print('Error: the result is not a list')
    else:
        print('[', end='')
        for i in range(len(y_f)-1):
            print("%.4f, "%y_f[i], end='')
        print("%.4f]"%y_f[len(y_f)-1])
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self, y): # y为当前状态，列表[S,I,R]
        S = y[0]
        I = y[1]
        R = y[2]
        return np.array([-self.beta * S * I, self.beta * S * I - self.gamma * I, self.gamma * I])
   
    def solve_with_control(self, x):
        # ------------------begin-----------------------
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1] - x[k]
            yp=y[k]+h*self.f(y[k])
            yc=y[k]+h*self.f(yp)
            y[k+1] = 1/2*(yp+yc)
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res






        # -------------------end------------------------

N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗） 
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常 [S0, I0, R0]
beta  = 1.0/N            # 平均传染率

for T in [5, 10, 20]:
    t   = range(0, T, 1) # 模拟T天的发展情况，单位时间为1天
    print("—————————————— 无出行控制下感染情况：——————————————")
    simulation = SIR(beta=beta, gamma=gamma, y0=y0)
    y_f = simulation.solve_with_control(t)
    myprint_list(y_f)

    for state in [2 , 5 , 10]:   # 出行控制的程度，对应感染率分别下降为原来的 1/state
        simulation = SIR(beta=beta*(1/state), gamma=gamma, y0=y0)
        y_f = simulation.solve_with_control(t)
        print("———出行控制后，传染率为原来的 1/{} 时,对应的感染情况：———".format(state))
        myprint_list(y_f)

第6关：绘制疫情发展趋势对比分析图

任务描述
每天的疫情数据代表着疫情的态势。通过对疫情数据可视化可以直观地反映疫情的态势。为有效防止疫情进一步扩散和蔓延，各地纷纷出台各种举措加强疫情防控工作。本关将在第 4 关（引入隔离机制）的基础上，绘制未隔离与隔离机制下新冠病毒发展趋势对比分析图。通过图示，可以直观看出为什么要引入隔离机制进行疫情防控？为什么要排查疑似？

测试说明
平台将运行用户补全的代码文件，并将存储的 src/step3/student/result.png 图像与标准答案图像比较，然后判断用户编写代码是否正确:

若画图正确，测试集将输出：祝贺！图片与预期输出一致;
否则，测试集将输出：图片与预期输出不一致，请继续努力！

在这里插入图片描述

#-*- coding : utf-8 -*-
# coding: utf-8
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

# 采用matplotlib作图时默认设置下是无法显示中文的，凡是汉字都会显示成方块。
# 实际上，matplotlib是支持unicode编码的，不能正常显示汉字主要是没有找到合适的中文字体。
from pylab import mpl
mpl.rcParams['font.sans-serif']= ['SimHei']
# 解决负号显示问题
matplotlib.rcParams['axes.unicode_minus']=False
        
class SIR():
    def __init__(self, beta, gamma, y0):
        self.beta = beta; self.gamma = gamma; self.y0 = y0  # 参数属性

    def f(self, y): #  y为当前状态，列表[S,I,R]
        S = y[0]
        I = y[1]
        R = y[2]
        return np.array([-self.beta * S * I, self.beta * S * I - self.gamma * I, self.gamma * I])
    
    def solve_with_quarantine(self, x):
        # ------------------begin-----------------------
        y = [self.y0] * len(x)
        for k in range(len(x)-1):
            h = x[k+1] - x[k]
            yp=y[k]+h*self.f(y[k])
            yc=y[k]+h*self.f(yp)
            y[k+1] = 1/2*(yp+yc)
        res=[]
        for i in range(len(y)):
            res.append(y[i][1])
        return res







        # -------------------end------------------------
        
N     = 1e8             # 武汉总人数：1000万人
gamma = 1/25            # 假设肺炎平均25天治愈（15天潜伏+10天治疗）
y0    = [N-1, 1, 0]     # 初始发病1人，其他人员正常 [S0, I0, R0]
beta = 1.0/N            # 平均传染率
t   = range(0, 60, 1)   # 模拟60天的发展情况，单位时间为1天
simulation = SIR(beta=beta, gamma=gamma, y0=y0)

# ------------------begin-----------------------
#1. 设置图形大小
plt.figure(figsize=(10, 8))

#2. 绘制曲线：横轴是时间/天，纵轴是感染人数
for gamma in [1/25,1/15,1/3]:
    simulation = SIR(beta=beta, gamma=gamma, y0=y0)
    y_f = simulation.solve_with_quarantine(t)
    plt.plot(t,y_f)


 
 

#3. 设置图题'未隔离与隔离机制下新冠病毒发展趋势对比分析'、
#   横轴标签'时间/天'、纵轴标签'人数'、
#   图列说明('未实施隔离', '确诊隔离', '疑似隔离', 分别对应三条曲线)
plt.title('未隔离与隔离机制下新冠病毒发展趋势对比分析')
plt.xlabel('时间/天')
plt.ylabel('人数')
plt.legend(['未实施隔离', '确诊隔离', '疑似隔离'])


# ------------------end-----------------------

# 设置y轴刻度
Vy = [1.0e6, 5.0e6] + [i * 1.e7 for i in range(11)]
plt.yticks(Vy, ['%d'%e for e in Vy])

plt.savefig('src/step6/student/result.png')
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值