dS/dt= -beta * S * I / N
dE/dt= beta * S * I / N - sigma * E 
dI/dt= sigma * E - gamma * I 
dR/dt= gamma * I

install.packages("deSolve")
library(deSolve)
N <- 1400000000 # 中国总人口数量
beta <- 0.13 # 病毒传播速率
gamma <- 0.995 # 治愈率（包括确诊与无症状感染者）
sigma <- 1/14 # 平均潜伏期
I <- 41699 #初始感染人数

library(deSolve)
# 定义 SEIRS 模型
seirs <- function(time, state, parms) {
  with(as.list(c(state, parms)), {
    # 模型方程
    dS <-  -beta*S*I/(S+E+I+R)
    dE <- beta*S*I/(S+E+I+R) - rho*E
    dI <- rho*E - gamma*I
    dR <- gamma*I
    # 返回微分方程组
    return(list(c(dS, dE, dI, dR)))
  })
}
# 设置模型参数
parms <- c(S = 1.4e9, E = 1, I = 41699, R = 0, rho = 1/14, gamma = 0.13, beta = 0.25)
# 设置初始状态
init <- c(S = 1.4e9, E = 1, I = 41699, R = 0)
# 设置起始时间为今天
start_time <- as.Date("2022-12-23")
# 模拟未来 N 天的趋势
N <- 365
times <- seq(start_time, by = "day", length.out = N)
times <- times <- as.numeric(times)
out <- as.data.frame(ode(y = init, times = times, func = seirs, parms = parms))
# 可视化结果
library(ggplot2)
df <- data.frame(time=as.Date(out$time, origin = "1970-01-01")
, S=out[,2], E=out[,3], I=out[,4], R=out[,5])
# 美化图表及设置线条属性
theme_set(theme_minimal(base_size = 14)) # 设置主题
line_size <- 1.5 # 设置线的宽度
line_type <- c('solid', 'dotted', 'dashed', 'dotdash') # 设置不同种类线条的样式
# 绘制图表
ggplot(df, aes(x = time)) +
  geom_line(aes(y = S, color = "Susceptible"), size = line_size, linetype = line_type[1]) +
  geom_line(aes(y = E, color = "Exposed"), size = line_size, linetype = line_type[2]) +
  geom_line(aes(y = I, color = "Infected"), size = line_size, linetype = line_type[3]) +
  geom_line(aes(y = R, color = "Recovered"), size = line_size, linetype = line_type[4]) +
  scale_color_manual(values = c("blue", "orange", "red", "green")) + # 自定义颜色
  labs(title = "SEIR Model for COVID-19 in China",
       x = "Time (days)",
       y = "Number of individuals",
       color = "Cases") +
  theme(legend.position = "bottom", legend.title = element_blank(),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 12),
        plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
        panel.border = element_rect(colour = "black", fill = NA),
        panel.grid.major = element_line(colour = "gray", linetype = "dotted", size = 0.5),
        panel.grid.minor = element_line(colour = "gray", linetype = "dotted", size = 0.3))
# 输出峰值时间和数量
peak_time <- df$time[which.max(df$I)]
peak_count <- max(df$I)
cat(paste("Peak time: ", peak_time, " days"))
cat("\n")
cat(paste("Peak count: ", peak_count, " individuals"))

Peak time:  2023-09-22
Peak count:  69027262.4873488