seir模型数学建模python_R示例：SEIR模型代码

说明：本例对常见的SEIR模型作了改变，由于DDE模型始终存在缺陷，故仍采用ODE建模。

实现环境：R 3.6.3+deSolve+ggplot2。

#-------------------------------------

#基于微分方程组建模

#dS/dt=-b_1*S*I/N-b_2*S*E/N

#dE/dt=b_1*S*I/N+b_2*S*E/N-k*E

#dI/dt=k*E-z*I

#dH/dt=z*I-d*H-m*H

#dR/dt=m*H

#dD/dt=d*H

#变量

#S易感人群

#E暴露潜伏

#I发病

#H入院

#R治愈

#D病故

#参数

#b_1,b_2接触感染概率

#k潜伏到发病的转化率

#z发病入院比例

#m治愈率

#d感染期死亡率

#--------------------------------------

library(deSolve)

library(ggplot2)

#模型

seir

with(as.list(c(state, pars)), {

dS

dE

dI

dH

dR

dD

dN

list(c(dS, dE, dI, dH, dR, dD, dN))

})

}

#参数

N

I_0

H_0

E_0

R_0

D_0

S_0

init

S = S_0,

E = E_0,

I = I_0,

H = H_0,

R = R_0,

D = D_0,

N = N

)

time

pars

b_1 = 0.16,

#发病者接触感染概率0.1575

b_2 = 0.79,

#潜伏期接触感染概率0.7875

k = 0.2,

#潜伏到发病的转化率

z = 0.25,

#发病入院比例

m = 0.04,

#治愈率

d = 0.0025

#感染期死亡率

)

#触发事件,S降低,10t至100%,20t至100%,90t维持控制

eventdat

var = c("S", "S", "S"),

time = c(10, 20, 90),

value = c(1, 1, 1),

method = c("mult", "mult", "mult")

)

#微分计算

res.seir

as.data.frame(lsoda(

y = init,

times = time,

func = seir,

parms = pars,

events = list(data = eventdat)

))

#可视化

ggplot(res.seir) +

geom_line(aes(x = time, y = N, col = "人口"), size = 1) +

geom_line(aes(x = time, y = S, col = "易感"), size = 1) +

geom_line(aes(x = time, y = E, col = "潜伏"), size = 2) +

geom_line(aes(x = time, y = I, col = "发病"), size = 2) +

geom_line(aes(x = time, y = H, col = "入院"), size = 2) +

geom_line(aes(x = time, y = R, col = "治愈"), size = 2) +

geom_line(aes(x = time, y = D, col = "病故"), size = 2) +

scale_color_manual(

name = "图例",

values = c(

"人口" = "purple",

"易感" = "blue",

"潜伏" = "orange",

"发病" = "red",

"入院" = "darkred",

"治愈" = "green",

"病故" = "black"

)

) +

scale_y_continuous("") +

labs(caption = "Note:All Hypothesis Reflecting No Reality")

结果：

2020-7-15 12:37:40 上传

下载附件 (5.23 KB)