# 一种SEIR模型的R代码

SEIR是很经典的模型，现在随着城市之间交通的日益频繁，考虑人口城市交通工具扩散的SEIR模型似乎更加符合实际。微分方程模型如下：

\begin{eqnarray*}

\frac{\mathrm{d}\vec{S}_t}{\mathrm{d}t} & = & -\vec{S}_t \odot \vec{\lambda}_t \\

\frac{\mathrm{d}\vec{E}_t}{\mathrm{d}t} & = & \vec{S}_t \odot \vec{\lambda}_t - \alpha \vec{E}_t \\

\frac{\mathrm{d}\vec{I}_t}{\mathrm{d}t} & = & \alpha \vec{E}_t - \gamma \vec{I}_t \\

\frac{\mathrm{d}\vec{R}_t}{\mathrm{d}t} & = & \gamma \vec{I}_t

\end{eqnarray*}

where $\vec{S}_t = (S_{1t},\dots,S_{nt})^T$ and likewise for $\vec{E}_t$, $\vec{I}_t$, and $\vec{R}_t$.  Furthermore,

$$\vec{\lambda}_t = \beta \left(\vec{I}/\vec{N} + (K \cdot (\vec{I_t}/\vec{N}))/\vec{N}\right)$$

K=diag(187)  #%% 城市之间交通流矩阵

N<<-china_population  #%% 城市人口向量

K<<-K

alpha=1/4

beta=2

gammaa=1/6

I0W = 1

param=c(beta, gammaa)

func = function(t, y, parms) {

n <- length(N) # number of cities

alpha=1/4

ymat = matrix(y, nrow=n, byrow=FALSE)

S = ymat[,1]

E = ymat[,2]

I = ymat[,3]

R = ymat[,4]

beta = parms[1]

gamma = parms[2]

lambda = beta * (I/N + (t(K/N)%*%I/N))

dS = -lambda * S

dE = lambda * S - alpha * E

dI = alpha * E - gamma * I

dR = gamma * I

list(c(dS, dE, dI, dR))

}

max_t=22

n = length(N) # number of cities

t = 0:max_t

# Initial infectives in Wuhan

I0 = rep(0, length(N))

I0[151] = I0W

y = c(N-I0, rep(0, length(N)), I0, rep(0, length(N)))

sim = deSolve::ode(y=y, times=t, func=func, parms=param)

t=sim[,1]

S=sim[,seq(2,len=n)]

E=sim[,seq(2+n, len=n)]

I=sim[,seq(2+2*n, len=n)]

R=sim[,seq(2+3*n, len=n)]

simulate=list(t=t, S=S, E=E, I=I, R=R)

write.csv(I,"mydata_I.csv")

write.csv(S,"mydata_S.csv")

write.csv(E,"mydata_E.csv")

write.csv(R,"mydata_R.csv")

plot(1:23, I[,151],type="o",

main="I")

grid()

box()

https://blog.sciencenet.cn/blog-306734-1219685.html

## 相关博文

GMT+8, 2023-5-28 17:53